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Positron  emission  tomography  (PET)  is  an  imaging  modality  that  provides  bio- 
chemical information  about  the  human  body.  Among  the  sources  of  error  in  PET, 
attenuation  is  the  one  that  impacts  image  quality  the  most.  The  standard  method 
for  correcting  attenuation  consists  of  the  following  two  step  procedure.  First,  a 
blank  scan  and  transmission  scan  are  used  to  estimate  the  survival  probabilities. 
Then,  the  data  for  each  pair  of  detectors  is  divided  by  the  corresponding  survival 
probability  estimate.  This  corrected  data  and  a  reconstruction  algorithm,  such  as 
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the  filtered  back-projection  algorithm,  are  then  used  to  generate  an  estimate  of 
the  emission  intensity,  which  is  displayed  as  an  image. 

In  this  dissertation,  two  iterative  methods  are  proposed  for  correcting  atten- 
uation errors  in  PET.  The  first  method  is  a  maximum  likelihood  (ML)  method 
that  simultaneously  estimates  the  emission  intensity  and  attenuation  density.  The 
second  one  is  a  least  squares  (LS)  reprojection  method,  where  the  survival  prob- 
abilities are  estimated  by  reprojecting  a  LS  estimate  of  the  attenuation  density. 
The  emission  intensity  is  estimated  using  the  expectation  maximization  maximum 
likelihood  (EMML)  algorithm  with  a  probability  matrix  that  has  been  modified  by 
the  survival  probability  estimates.  Both  methods  are  guaranteed  to  produce  non- 
negative  estimates  of  the  emission  intensity  and  the  attenuation  density,  provided 
that  the  initial  estimates  are  nonnegative. 

From  simulation  studies,  it  appears  that  the  major  advantage  of  the  proposed 
methods  over  the  standard  method  with  the  filtered  back-projection  algorithm 
is  that  they  (1)  have  less  noise,  (2)  do  not  have  artifacts  outside  the  region  of 
interest,  and  (3)  they  provide  high  quality  images  even  for  short  transmission 
scans  (less  than  10  minutes). 
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CHAPTER  1 
INTRODUCTION 

For  both  clinical  and  research  applications,  medical  imaging  modalities  are 
used  to  obtain  important  information  on  the  functionality  and  anatomy  of  the 
human  body.  Imaging  modalities  such  as  X-ray  computed  tomography  and  mag- 
netic resonance  imaging  provide  detailed  images  of  the  anatomical  structure  of 
the  body.  However,  in  many  applications,  it  is  more  beneficial  to  have  a  modality 
that  determines  functional  information  such  as  blood  flow  rates,  oxygen  utiliza- 
tion rates,  and  glucose  utilization  rates  within  the  region-of-interest  (ROI)  (e.g., 
brain,  heart)  [28].  For  example,  glucose  utilization  rates  are  used  in  oncology  to 
diagnose  tumors  and  to  determine  the  efficacy  of  a  treatment. 

Functional  information  on  the  human  body  can  be  obtained  using  (1)  positron 
emission  tomography  (PET),  (2)  single-photon  emission  computed  tomography 
(SPECT)  [16],  and  (3)  functional  magnetic  resonance  imaging  [51].  Among  them, 
PET  is  probably  the  most  flexible  imaging  modality.  For  example,  PET  can 
image  which  parts  of  the  brain  are  functioning  while  a  subject  is  performing  a 
speaking  tasking.  This  is  not  possible  in  functional  magnetic  resonance  imaging 
because  the  scanners  are  too  loud  [55].  An  advantage  of  PET  over  SPECT  is 
that  the  radionuclides  in  PET  can  be  incorporated  into  molecules  (or  analogs)  of 
biological  interest  [56]. 


In  this  dissertation,  we  propose  two  methods  for  correcting  a  certain  error  in 
PET  called  attenuation.  The  outline  of  the  dissertation  is  as  follows.  An  overview 
of  PET  is  provided  in  the  remainder  of  Chapter  1.  In  Chapter  2,  we  discuss  the 
physics  behind  PET,  introduce  a  mathematical  model  for  the  data,  and  review 
existing  attenuation  correction  methods.  We  introduce  the  proposed  maximum 
likelihood  and  least-squares  reprojection  methods  in  Chapter  3  and  Chapter  4, 
respectively.  In  Chapter  5,  we  investigate  the  performance  of  these  methods  using 
synthetic  and  real  data.  Finally,  we  discuss  our  conclusions  and  future  research 
direction  in  Chapter  6. 


1.1      Overview:   Positron  Amission  Tomography 

PET  is  exceptional  in  its  capability  to  quantitatively  determine  functional 
information  within  a  ROI.  Since  George  Hevesy  won  a  Nobel  Prize  for  his  pio- 
neering research  on  radionuclides  in  1943,  numerous  researchers  have  investigated 
methods  to  obtain  functional  information  on  the  human  body.  However,  the  first 
PET  machine  was  not  built  until  1975  by  Michael  Phelps,  Edward  Hoffman,  and 
their  colleagues  [42].  Since  then,  scientists  and  engineers  have  made  significant 
advances  in  the  hardware,  radionuclides,  and  reconstruction  algorithms  used  for 
PET.  In  Fig.  1.1  is  a  PET  scanner  of  the  Quest  series  from  the  GE  Medical 
Systems. 

We  first  introduce  the  major  components  of  a  PET  scan,  which  are  as  follows: 

1.  A  compound  is  labeled  with  a  radionuclide. 


Figure  1.1.   Quest  250  from  GE  Medical  Systems  (courtesy  of  GE  Medical  Sys- 
tems) 


2.  The  labeled  compound  or  radiopharmaceutical  is  introduced  into  the  sub- 
ject's body. 

3.  The  data  are  collected. 

4.  An  image  is  produced  that  represents  the  concentration  of  the  radiophar- 
maceutical within  the  ROI. 

Radionuclides  in  PET  are  (1)  incorporated  into  compounds  that  are  naturally 
used  by  the  body,  (2)  emit  positrons,  and  (3)  have  short  half-lives.  A  radionuclide's 
half-life  must  be  sufficiently  long  so  that  there  is  enough  data  to  get  quality 
images  [38].  The  most  common  radionuclides  in  PET  include  nC,  150,  and 
13N,  which  are  used  to  label  glucose,  water  and  amino  acid,  respectively.  The 
most  widely  used  radiopharmaceutical  is  18F-labeled  deoxyglucose  (FDG),  which 
is  used  to  study  glucose  metabolism.  Because  of  the  short  half-lives  of  positron- 
emitting  radionuclides  (2  minutes,  10  minutes,  20  minutes,  and  110  minutes  for 
nC,  150,  13N,  and  18F,  respectively),  they  need  to  be  produced  on-site  using  a 
cyclotron.  The  cyclotron  creates  radionuclides  by  smashing  the  atomic  shell  of 
certain  molecules  with  protons.  Cyclotrons  are  expensive  and  costly  to  maintain 
because  they  require  at  least  two  operators. 

SPECT  uses  radionuclides,  such  as  99mTc,  201TI,  and  123I,  which  have  half- 
lives  that  range  from  10  hours  to  several  days.  Because  of  their  long  half-lives, 
radionuclides  in  SPECT  can  be  made  in  advance  and  shipped  to  hospitals  so  an 
on-site  cyclotron  is  unnecessary.  However,  long  half-lives  mean  that  patients  are 
exposed  to  radioactivity  for  long  periods. 


After  being  administered  to  a  subject,  the  labeled  compound  will  localize  in  the 
ROI.  For  example,  FDG  is  absorbed  preferentially  by  the  brain  because  glucose 
is  a  source  of  energy  [11].  Another  example  is  the  use  of  150  based  radiopharma- 
ceuticals to  index  injured  tissues  by  monitoring  the  oxygen  metabolism  of  tissues 
[52]. 

1.2     Emission  Data 

As  mentioned  previously,  the  radiopharmaceutical  localizes  in  the  ROI.  When 
the  radionuclide  decays,  it  emits  positrons  which  travel  a  short  distance  (3-4  mm) 
before  annihilating  with  nearby  electrons.  The  areas  with  a  high  concentration  of 
the  radiopharmaceutical  emit  more  photons  than  the  areas  of  low  concentration. 
When  a  positron  annihilates  with  an  electron,  two  high-energy  (511  keV)  photons 
are  produced  that  fly  off  at  an  almost  180-degree  angle  from  each  other.  This 
annihilation  process  is  shown  in  Fig.  1.2.  The  photons  go  through  the  subject  and 
are  registered  by  rings  of  detectors  surrounding  the  subject.  When  two  detectors 
detect  a  pair  of  photons  within  a  short  period  of  time,  these  photons  are  said  to 
be  in  coincidence  (see  Fig.  1.2).  For  each  pair  of  detectors,  there  is  a  count  that 
corresponds  to  the  number  of  observed  coincidences.  The  emission  data  is  the 
collection  of  these  counts. 

From  the  emission  data,  we  can  estimate  the  concentration  of  the  radiophar- 
maceutical using  a  reconstruction  algorithm.  A  computer  is  then  used  to  display 
images  where  the  intensity  is  proportional  to  the  estimated  concentration  of  the 
radiopharmaceutical.  In  the  next  section,  some  clinical  and  research  applications 
of  PET  will  be  presented. 


Figure  1.2.  A  positron  travels  a  short  distance  and  annihilates  with  a  nearby 
electron  to  produce  a  pair  of  photons  that  fly  off  at  180-degree  angle.  These 
photons  are  detected  by  a  ring  of  detectors 


1.3     Clinical  and  Research  Applications 

PET  is  widely  used  in  pathology  and  oncology  because  it  can  quantitatively 
measure  blood  flow,  and  utilization  rates  of  glucose,  CO2,  and  fatty  acid  ammo- 
nia. PET  provides  functional  information  about  tissues  or  organs,  which  cannot 
be  obtained  by  modalities  such  as  magnetic  resonance  imaging  and  computed  to- 
mography [52].  It  is  also  a  powerful  tool  for  investigating  brain  functions  such 
as  emotion,  thinking,  vision,  and  speech  [5,  30].  Since  the  major  energy  source 
of  the  brain  is  glucose,  FDG  is  used  to  map  the  brain's  energy  consumption.  In 
other  words,  if  a  specific  region  of  the  brain  has  increased  functionality,  then  the 
regional  blood  flow  increases  to  provide  the  necessary  energy.  As  a  result,  the  con- 
centration of  the  radionuclide  used  in  FDG  increases  in  that  region.  Using  PET 
to  map  the  concentration  of  the  radionuclide,  we  can  observe  the  higher  function- 
ality in  that  region  of  the  brain.  Based  on  the  concepts  discussed  above,  scientists 
have  used  PET  images  to  assess  the  brains  of  normal  children  and  children  under 
abuse.  They  found  that  the  children  who  were  abused  had  slower  developments 
in  brain  functions  such  as  speech  and  vision  as  compared  to  normal  children  [7]. 

PET  has  clinical  applications  as  well.  For  example,  PET  has  the  ability  to 
measure  the  size  of  infarcts  in  the  heart  after  a  coronary  attack,  to  diagnose  dis- 
eases such  as  Alzheimer's  [45]  and  epileptic  seizures  [36],  and  to  survey  the  effects 
of  treatments  such  as  radiation  therapy,  chemotherapy,  and  hormone  therapy.  In 
addition,  whole-body  PET  is  used  to  spot  metastases,  which  is  very  difficult  for 
other  modalities  because  of  the  small  sizes  and  the  random  occurrences  of  metas- 
tases inside  cancer  patients'  bodies.    The  effect  of  the  chemotherapy  on  cancer 
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patients  can  be  evaluated  by  visualizing  the  sizes  and  locations  of  metastases  on 
PET  images  [1]. 

PET  is  also  used  in  psychological  studies.  For  instance,  it  was  used  in  the 
early  1990s  several  times  by  courts  in  California  to  determine  whether  neural 
malfunctions  in  the  brain  caused  a  suspect's  violent  behavior.  Another  interesting 
experiment  was  the  monitoring  of  six  bank  officials'  brain  functions  when  they 
were  exposed  to  a  video  of  an  armed  bank  robbery.  This  was  done  to  find  out  the 
connection  between  stress  and  visual  re-experience  [22]. 


1.4     Advantages  and  Disadvantages 

PET  is  a  better  imaging  modality  than  SPECT  in  the  sense  of  offering  better 
efficiency,  sensitivity,  and  accuracy  as  compared  to  SPECT  [11].  Also,  unlike 
SPECT,  the  radionuclides  used  in  PET  are  easily  incorporated  into  compounds, 
such  as  glucose,  that  are  naturally  used  by  the  body.  Using  different  radionuclides, 
PET  is  able  to  image  different  biochemical  activities  while  functional  magnetic 
resonance  imaging  only  can  image  hydrogen  activity.  Additionally,  PET  can  track 
chemical  changes  that  do  not  affect  blood  flow  enough  for  functional  magnetic 
resonance  imaging  to  measure. 

Unfortunately,  PET  has  some  major  disadvantages  such  as  its  high  cost,  high 
complexity,  and  the  requirement  of  highly  skilled  personnel.  The  cost  of  a  PET 
scanner  is  similar  to  a  magnetic  resonance  imaging  scanner.  However,  the  cy- 
clotron used  to  generate  the  radionuclides  is  also  very  expensive.  Additionally, 
at  least  three  skilled  professionals  -  one  to  maintain  the  scanner,  one  to  run  the 


cyclotron  and  another  one  to  produce  radionuclides  -  are  needed  to  operate  a 
PET  laboratory.  For  SPECT,  it  is  not  necessary  to  create  the  radionuclides  on 
site  because  the  radionuclides  used  in  SPECT  can  be  stored.  Even  though  its 
image  quality  is  not  as  good  as  that  of  PET,  SPECT  is  more  widely  available  in 
practice  because  it  is  significantly  cheaper. 


CHAPTER  2 
BACKGROUND 

In  this  chapter,  we  review  the  basic  physics  of  PET,  introduce  the  Poisson 
model,  and  discuss  sources  of  errors  in  the  data.  Since  the  focus  of  this  dissertation 
is  attenuation  correction,  we  also  review  existing  attenuation  correction  methods. 


2.1     Basic  Physics  of  PET 

As  discussed  in  Chapter  1,  the  subject  is  given  a  radiopharmaceutical  that 
emits  positrons.  When  a  positron  annihilates  with  a  nearby  electron,  two  photons 
are  produced  that  depart  at  an  angle  near  180  degrees.  The  region  of  interest 
(ROI)  is  surrounded  by  rings  of  detectors  that  detect  the  photons.  We  will  use 
the  term  tube  to  refer  to  a  pair  of  detectors.  When  a  pair  of  photons  is  incident 
on  a  tube  within  a  small  time  interval,  a  counter  is  increased  (see  Fig.  2.1).  The 
emission  data  are  the  counts  associated  with  the  various  tubes. 

The  detectors  are  shielded  from  radiation  coming  from  outside  the  region  of 
interest  by  heavy  lead  rings  at  the  front  and  back  of  the  gantry.  In  order  to 
switch  between  the  slice-collimated  mode  and  the  three-dimensional  mode,  most 
commercial  PET  scanners  contain  spherical  thin  interplane  rings  called  septa. 
When  the  scanner  is  in  slice-collimated  mode,  the  septa  only  allows  in-plane  tubes 
to  detect  photons.   When  a  scanner  operates  in  full  three-dimensional  mode,  the 
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Ring  of  Detectors 


Figure  2.1.  Annihilation  Process 
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septa  are  withdrawn  to  permit  photons  to  be  registered  by  in-plane  tubes  and 
cross-plane  tubes  [37]. 

A  detector  module  consists  of  an  array  of  crystals  coupled  to  several  photomul- 
tiplier  tubes.  In  Fig.  2.2,  a  detector  module  of  the  Scanditronix  PC4096-15WB 
whole-body  PET  scanner  is  shown.  A  detector  ring  contains  several  hundred  de- 
tector modules  (e.g.,  384  for  the  Siemens-CTI  ECAT  EXACT  921  scanner).  When 
a  crystal  is  hit  by  a  511  keV  photon,  it  interacts  with  the  photon  and  generates  nu- 
merous light  photons.  Then,  the  light  guide  connected  to  the  crystal  transfers  the 
light  photons  to  the  nearest  photomultiplier  tube  (PMT).  The  PMTs  collect  light 
photons  and  convert  them  into  electrical  signals  that  are  stored  in  a  computer. 
The  outputs  of  the  PMTs  are  decoded  to  determine  which  detector  module  was 
hit  by  a  photon.  The  most  widely  used  crystal  material  is  Bismuth-Germanate 
because  it  has  good  properties  when  it  interacts  with  511  keV  photons  such  as  a 
high  light  photon  producing  rate,  short  decay  time,  and  a  fast  energy  rise  time 
[26]. 

Once  the  data  is  collected,  algorithms  are  used  to  reconstruct  the  PET  im- 
ages. The  most  popular  method  is  the  filtered  back-projection  algorithm  [45].  It 
is  extremely  fast  and  easy  to  implement,  but  the  images  are  noisy,  have  nega- 
tive values,  and  contain  artifacts  outside  the  ROI.  Alternative  approaches  include 
statistically  based  methods  such  as  maximum  likelihood  [28,  46],  maximum  a 
posterior  [20,  21,  30,  31],  least  squares  [15,  40,  24,  46,  49,  50],  and  weighted  least- 
squares  [3,  12,  17]  methods.  An  advantage  of  statistically  based  methods  is  that 
they  account  for  the  errors  in  the  data  and  do  not  have  artifacts  outside  the  ROI. 
A  disadvantage  of  them  is  that  they  are  computationally  demanding. 
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Light  Guide 


Crystals 


Figure  2.2.     Scanditronix  PC4096-15WB  detector  module:    16  crystals,  1  light 
guide,  and  2  dual  cathode  PMTs 
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2.2     Poisson  Model 

Since  the  data  collection  in  PET  is  actually  a  photon  counting  process,  it  is 
reasonable  to  assume  that  the  emission  process  is  a  Poisson  process  with  unknown 
mean  [38].  Based  on  this  assumption,  Shepp  and  Vardi  proposed  a  model  called 
the  Poisson  model  [46]  in  1982.  Their  model  assumes  that  the  data  is  error  free 
and  obeys  a  Poisson  process. 

In  the  Poisson  model,  the  ROI  is  discretized  into  a  grid  of  J  voxels.  Let  Nj 
denote  the  number  of  emissions  in  voxel  j,  where  Nj  is  assumed  to  be  Poisson  with 
unknown  mean  Xj  =  E[Nj].  Let  Y{  denote  the  number  of  coincidences  detected  in 
the  ith  tube,  where  Y{  is  assumed  to  be  Poisson  with  mean  d,  =  E[Yi\  and  /  is  the 
number  of  tubes.  It  is  also  assumed  that  (1)  {Nj}{  are  independent,  (2)  {Yj}{ 
are  independent,  and  (3)  the  probability  that  an  emission  in  voxel  j  is  detected  in 
the  ith  tube,  P,j,  is  known  for  all  i,  j.  Under  these  assumptions,  it  can  be  shown 
that  the  mean  of  the  ith  tube  is  given  by 

j 

dr     =    Y,X,Pij  (2.1) 

.1  =  1 

=    {Px)i,       t  =  l,2,... ,/,  (2.2) 

where  the  (i,j)th  element  of  P  is  PtJ  and  (Pa;),  denotes  the  ith  element  of  the 
vector  Px. 

Given  the  data,  3/1,2/2,  ■  ■  ■  ,Vi,  where  ?/,  is  an  observation  of  Y{,  the  goal  is  to 
estimate  x  =  [x\  x?  ■  ■  •  Xj]  .  In  this  way,  an  estimate  for  the  number  of  photons 
emitted  in  each  voxel  is  obtained.  Throughout  the  remainder  of  this  dissertation, 
x  will  be  referred  to  as  the  emission  intensity. 
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Shepp  and  Vardi  used  the  expectation  maximization  algorithm  [16]  to  compute 
a  ML  estimate  of  the  emission  intensity.  This  algorithm,  which  is  referred  to  as 
the  expectation  maximization  maximum  likelihood  algorithm,  is  briefly  discussed 
in  the  Appendix. 

In  the  next  section,  we  discuss  sources  of  errors  in  PET. 


2.3     Errors  in  Positron  Emission  Tomography 

In  reality,  PET  data  are  corrupted  by  errors.  In  this  section,  we  go  through 
most  of  the  important  sources  of  error.  Of  all  the  errors,  the  one  that  has  the 
most  effect  on  image  quality  is  attenuation  [23].  Since  we  focus  on  developing  new 
algorithms  to  correct  attenuation,  we  provide  an  in  depth  discussion  of  attenua- 
tion. 


2.3.1     Background  on  Attenuation 

Attenuation  is  the  main  error  source  in  PET  emission  image  reconstruction 
[34].  The  511  keV  photons  produced  by  the  annihilations  are  attenuated  by  tissue, 
bone,  and  air.  Different  mediums  have  different  attenuation  effects.  Before  going 
into  the  details  of  attenuation  correction  for  PET,  we  introduce  the  physics  behind 
attenuation. 

The  main  sources  of  attenuation  in  PET  are  absorption  and  Compton  scatter- 
ing.  Photons  are  rarely  absorbed  because  they  carry  511  keV  of  energy.   Hence, 
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Figure  2.3.  Attenuation  due  to  absorption 
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attenuation  due  to  absorption  is  negligible.  Fig.  2.3  is  a  depiction  of  photon  ab- 
sorption. Compton  scattering  occurs  when  the  path  of  the  photon  is  deflected  by 
an  outer  shell  electron  [37].  Photons  in  PET  are  more  likely  to  undergo  Compton 
scattering  than  absorption.  Most  of  the  photons  that  undergo  Compton  scatter- 
ing are  not  detected.  Fig.  2.4  shows  a  photon  undergoing  Compton  scattering. 
The  total  effect  of  photon  absorption  and  Compton  scattering  is  referred  to  as 
attenuation. 

The  survival  probability  along  a  line  L,  a,  is  the  probability  that  a  photon  is 
not  attenuated  along  that  line.  It  can  be  shown  that  the  survival  probability  is 
equal  to 

a  =  e-ft'uU,  (2.3) 

where  fi  is  the  attenuation  density  [29]. 

Since  eq.(2.3)  is  a  fundamental  equation  in  attenuation,  we  provide  an  expla- 
nation of  it.  Consider  Fig.  2.5,  which  shows  a  line  source  emitting  photons  toward 
an  object  with  attenuation  density  \i.  We  assume  that  N  photons  incident  on  a 
point  within  the  object  travel  a  short  distance,  A/,  and  only  N  +  AN  photons 
survive.  Obviously,  AN  is  a  negative  number.  We  may  use  the  following  equation 
to  describe  the  above  process: 


fir- 
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Figure  2.4.  Attenuation  due  to  scattering 
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Figure  2.5.  Basic  concept  of  attenuation 
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As  A/  and  AN  go  to  zero,  we  get  the  differential  equation 

&~" 

which  can  be  solved  by  integrating  along  the  line  L.  Doing  so,  we  have  that 

/       ^f  =  -  / jull ,  (2.6) 

JNtn      N  Jl 

where  7V,-„  is  the  number  of  photons  that  go  into  the  object  and  Nout  is  the  number 
of  photons  that  come  out  of  the  object.  From  the  above  equation,  we  get 

\nNout-\nNin  =  -J  fidl.  (2.7) 

It  follows  from  eq.  (2.7)  that 


a; 


out 


srC"/*"*.  (2.8) 


Nin 

Finally,  eq.  (2.3)  follows  by  noting  that  the  left  side  of  (2.8)  is  the  probability 
that  a  photon  along  the  line  L  is  not  attenuated.  In  the  next  several  sections,  we 
discuss  some  other  errors  in  PET. 


2.3.2     Accidental  Coincidences 

If  two  photons  from  different  annihilations  are  registered  as  a  coincidence  by 
a  tube,  then  an  error  occurs.  This  kind  of  error  is  referred  to  as  an  accidental 
coincidence  [22].  Fig.  2.6,  shows  an  example  of  an  accident  coincidence.  Two 
annihilations  happen  inside  the  ROI.  However,  in  both  annihilations,  a  photon 
flies  out  of  the  detector  ring.  The  remaining  photons  are  detected  by  a  tube 
within  the  coincidence  window  and  incorrectly  recorded  as  a  coincidence. 
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Out  of  plane 
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Figure  2.6.  Accidental  coincidence  in  PET 
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The  accidental  coincidence  rate  is  proportional  to  the  concentration  of  the 
radionuclide  in  the  ROI.  Usually,  the  accidental  coincidences  are  corrected  by 
subtracting  an  estimate  of  them  from  the  observed  photon  counts.  Accidental 
coincidences  are  estimated  either  by  using  an  estimate  of  the  single  events  to 
compute  them  or  using  a  delayed  coincidence  circuitry  to  measure  them  [22]. 
This  subtraction  method  may  produce  negative  counts  in  the  corrected  data, 
which  means  that  it  is  no  longer  Poisson.  Consequently,  several  methods  have 
been  proposed  to  resolve  this  problem  [2,  28,  41].  A  discussion  can  be  found  in 
Anderson  et  al.  [4]  and  Ollinger  et  al.  [37].  In  the  next  section,  another  error 
source,  scatter,  is  reviewed. 


2.3.3     Scatter 

In  a  scattered  event,  one  or  both  of  the  511  keV  photons  created  in  an  an- 
nihilation are  scattered  but  they  still  intersect  a  detector  ring.  As  a  result,  the 
annihilation  event  is  registered  by  the  wrong  tube.  An  example  of  an  error  due 
to  scatter  is  shown  in  Fig.  2.7. 

Since  photons  lose  energy  after  they  are  scattered,  it  is  very  easy  to  distinguish 
them  from  unscattered  photons  by  measuring  the  energy  they  pass  to  the  detector 
crystals  [8].  This  technique  can  be  used  to  correct  the  scattered  events  by  setting 
a  threshold  to  reject  photons  with  energy  sufficiently  less  than  511  keV.  Scatter 
between  different  slices  in  whole-body  PET  is  minimized  by  the  septa  (see  Section 
2.1). 
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Figure  2.7.  Scatter  in  PET 
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2.3.4     Other  Errors 

In  reality,  positrons  travel  a  short  distance  (3-4  mm)  before  annihilating  with 
an  electron  and  photons  shoot  off  at  an  angle  less  than  180°.  Also,  photons  can 
penetrate  the  crystal  they  hit  and  leak  energy  to  a  neighboring  tube  [26].  Another 
type  of  error,  called  detector  inefficiency,  occurs  when  a  detector  fails  to  detect 
an  incident  photon  [17].  Fortunately,  most  of  the  time,  these  errors  are  small  and 
usually  can  be  ignored  or  corrected. 

In  order  to  correct  for  attenuation,  certain  scans  called  blank  scans  and  trans- 
mission scans  are  used  to  estimate  the  survival  probability  along  each  tube.  In 
the  next  section,  we  provide  a  brief  discussion  of  blank  scans  and  transmission 
scans. 


2.4     Transmission  Tomography  in  PET 

In  order  to  correct  for  attenuation,  a  blank  scan  and  transmission  scan  are 
performed.  For  the  blank  scan,  a  scan  is  done  without  a  subject.  The  external 
photon  source,  as  seen  in  Fig.  2.8,  rotates  along  the  interior  of  the  detector  ring. 
Let  6,  denote  the  number  of  observed  coincidences  in  tube  i  during  the  blank  scan, 
i  =  1,2,...,/.  The  blank  scan  lasts  about  an  hour  to  more  than  2  hours.  It  is 
performed  at  the  beginning  of  every  day. 

After  positioning  the  subject  in  the  scanner,  a  scan  is  performed,  which  is 
referred  to  as  a  transmission  scan.  Like  the  blank  scan,  the  transmission  scan  uses 
a  rotating  external  source  of  photons  (see  Fig.  2.9).  The  photon  counts  collected 
in  a  transmission  scan  are  the  ones  that  survive  the  attenuation  effect  of  the  ROI. 
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Figure  2.8.  Blank  scan 
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Let  ri  denote  the  number  of  observed  coincidences  during  the  transmission  scan 
for  tube  i,  i  =  1, 2, . . . ,  /.  A  transmission  scan  needs  to  be  performed  each  time 
a  new  patient  is  put  in  the  PET  scanner.  The  duration  of  a  transmission  scan 
ranges  from  about  2  to  25  minutes.  For  the  safety  and  comfort  of  the  patient,  the 
transmission  scans  are  kept  as  short  as  possible,  (patient  should  not  be  exposed 
to  the  radioactive  source  too  long). 

Intuitively,  it  is  reasonable  to  estimate  the  survival  probability  by  taking  the 
ratio  of  the  transmission  scan  count  per  unit  time  and  the  blank  scan  count  per 
unit  time.  Hence,  an  estimate  of  the  survival  probability  in  tube  i,  a,-,  is 

a,  =  4^  .  (2-9) 

where  rj  and  rr  are  the  durations  of  the  blank  scan  and  the  transmission  scan, 
respectively.  The  data  can  then  be  corrected  for  attenuation  by  dividing  the  data 
by  the  estimated  survival  probabilities  (i.e.,  ^).  This  approach  is  simple  and 
intuitively  appealing,  but  it  does  not  work  well  if  the  transmission  scan  is  too 
short.  One  way  to  resolve  that  problem  is  to  modify  the  Poisson  model  instead  of 
directly  correcting  the  data.  In  next  section,  we  introduce  the  modified  Poisson 
model. 


2.5     Modified  Poisson  Model 

In  this  section,  we  discuss  an  extension  of  the  Poisson  model  for  the  case 
of  attenuation  errors  only  [28].  Before  discussing  this  model,  we  introduce  the 
concept  of  thinning  a  Poisson  process. 
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Figure  2.9.  Transmission  scan 


28 


Result  1:  Let  N  be  a  Poisson  random  variable  with  mean  x.  Suppose  a  random 
variable  X  is  created  by  thinning  N  with  probability  p.  Then  X  is  Poisson 
with  mean  px. 

Proof: 


P(X  =  k)    =     J2  P(X  =  k\N  =  m)P(N  =  m) 

m  =  k 
oo 

=     E(?)Pk(l-P)m-kp(N  =  m) 


i=k 

^  m! 


^kk\(m-k)\P[        V)  m! 

e~y    ~    (l-p)m-kxm 


Replacing  m  —  k  with  n,  yields 


P(X  =  t)  .  e^gd-pr^ 


&!     „^0  n! 


e-^(xp)fc   -   ((l-p)s)" 


_     e-'(xp)k  (1_p)x 
k\ 
e~xp(xp)k 
~  k\         ' 

which  shows  that  X  is  a  Poisson  random  variable  with  mean  px. 

We  can  apply  Result  1  to  model  the  effect  of  attenuation  in  PET  (see  also 
[41]).  Let  pj  be  the  attenuation  density  within  the  jih  voxel,  where  it  is  assumed 
that  the  attenuation  density  within  each  voxel  is  constant.  Let  L  be  an  /  x  J 
matrix  whose  (i,j)th  element,  L,j,  is  the  intersection  length  of  the  ith  tube  and 
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the  jih  voxel.  Then,  it  follows  from  eq.(2.3)  that  the  probability  that  a  photon 
pair  is  detected  by  tube  i  is  approximately  equal  to 

at  ±  e-lLM'  ,  (2.10) 

where  (Lfi)i  —  Yli=i  LijPji  i  =  1,2, ...,/.  Because  of  attenuation,  the  data  in 
tube  i,  Yi,  is  created  by  thinning  the  error  free  data  in  tube  i  by  the  survival 
probability  a,.  By  Result  1,  it  follows  that  Y{  is  Poisson  with  unknown  mean 

E[Yi]  =  ai{Px)i  .  (2.11) 

With  this  modified  Poisson  model,  several  statistical  reconstruction  algorithms 
have  been  proposed. 


2.6     Existing  Attenuation  Correction  Methods 

In  this  section,  we  review  existing  attenuation  correction  methods  including 
those  based  on  the  modified  Poisson  model. 


2.6.1     Standard  Correction  Method 

The  standard  way  to  compensate  for  the  attenuation  is  to  divide  the  emission 
data  by  the  estimate  of  the  survival  probabilities.  Specifically,  the  corrected  data 
for  the  iih  tube,  y,,  is  given  by 

fc=^-,i  =  l,2,...,/,  (2.12) 
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where  d;  is  given  by  eq.  (2.9). 

This  method  is  easy  to  implement  and  works  fairly  well  for  long  transmission 
scans  (20-25  minutes),  for  short  transmission  scans  (less  than  10  minutes),  it 
produces  noisy  emission  images.  The  conventional  way  this  problem  is  addressed 
is  to  smooth  the  transmission  sinogram.  Smoothing  reduces  the  noise,  but  it  also 
lessens  the  quantitative  accuracy  of  the  images  [34]. 


2.6.2     Reprojection  Methods 

Given  the  blank  scan  data  and  the  transmission  scan  data,  an  estimate  of  the 
attenuation  density  is  obtained  [36,  41].  Then,  the  estimate,  ji,  is  "reprojected" 
to  obtain  an  estimate  of  the  survival  probabilities: 

^-e-Eii^fo   ,i  =  1,2,...,/.  (2.13) 

The  estimated  survival  probabilities  are  used  to  correct  the  data  or  modify  the 
probability  matrix,  P.  In  addition  to  improving  the  accuracy  of  the  estimated 
survival  probabilities,  the  transmission  images  provide  information  that  can  be 
used  to  position  patients,  estimate  scattering  errors,  and  to  provide  anatomical 
landmarks  that  cannot  be  seen  in  emission  images  [36]. 

The  reprojection  method  provides  a  better  estimate  of  the  survival  probabil- 
ities. However,  the  counting  noise  from  the  transmission  image  introduces  addi- 
tional error  into  the  emission  reconstruction.  In  regions  such  as  the  thorax,  where 
the  attenuation  medium  is  heterogeneous,  this  problem  cannot  be  resolved  by 
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simply  smoothing  the  estimated  survival  probabilities  because  it  introduces  bias 
[13].  Since  reprojection  methods  estimate  the  attenuation  density,  they  are  more 
computationally  demanding  than  the  standard  correction  method.  In  the  next 
section,  we  discuss  methods  that  simultaneously  estimate  the  attenuation  density 
for  survival  probabilities  and  the  emission  intensity. 


2.6.3     Simultaneous  Estimation  Methods 

In  [41],  Politte  and  Snyder  proposed  constrained  maximum  likelihood  recon- 
struction algorithms  by  modifying  the  probability  matrix  under  the  assumption  of 
known  survival  probabilities.  Their  methods  are  not  practical  since  the  survival 
probabilities  must  be  estimated.  The  transmission  data  provides  noisy  estimate 
of  the  survival  probabilities  because  the  count  is  pretty  low  in  real  applications. 
Since  the  emission  data  are  more  reliable  than  the  transmission  data,  in  the  sense 
of  having  higher  data  counts,  simultaneous  estimation  methods  may  provide  a 
more  accurate  estimate  of  the  emission  intensity. 

In  [11],  a  simultaneous  estimation  method  was  proposed  to  estimate  the  atten- 
uation density  and  the  emission  intensities  for  SPECT  using  the  Cyclic  Subgradi- 
ent  Projection  (CSP)  algorithm.  After  the  Poisson  model  was  proposed,  Fessler, 
Clint  home  and  Rogers  proposed  a  joint  maximum  likelihood  method  that  esti- 
mates the  emission  intensity  and  survival  probabilities  [13,  18].  Lu  proposed  a 
MAP  algorithm  to  estimate  the  emission  intensity  and  the  attenuation  density  for 
SPECT  [32].  However,  prior  information  about  the  attenuation  medium  must  be 
known. 
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To  simultaneously  estimate  the  attenuation  density  and  the  emission  intensity, 
we  propose  a  new  method.  In  the  next  chapter,  we  discuss  this  proposed  method. 


CHAPTER  3 
MAXIMUM  LIKELIHOOD  METHOD 

In  this  chapter,  we  present  an  attenuation  correction  method  that  jointly  esti- 
mates the  attenuation  density  and  emission  intensity  using  a  maximum  likelihood 
(ML)  method.  Other  ML  methods  have  been  proposed  [13,  18],  but  they  have 
the  disadvantage  of  requiring  more  unknowns  than  data  points.  By  contrast, 
the  proposed  method  has  fewer  unknowns  than  data  points  for  typical  scanners. 
Specifically,  it  estimates  2J  unknowns  from  /  data  points,  where  J  and  /  are  the 
number  of  the  voxels  and  the  number  of  the  tubes,  respectively.  A  common  choice 
for  the  number  of  voxels  is  J  =  16384  (i.e.,  J  =  128  x  128)  and,  for  the  Siemens- 
CTI  ECAT  EXACT  921  PET  scanner,  the  number  of  the  tubes  is  /  =  73536. 


3.1     Maximum  Likelihood  Estimator 

As  discussed  in  Section  2.2,  the  emission  data  for  the  ith  tube,  i/,-,  is  assumed 
to  be  an  observation  of  a  random  variable,  Y^,  where  Yi  is  Poisson  with  mean 
{Px)ie-(LM<  (recall,  (Px){  =  EJj=ixJPiJ  and  i1^  =  Ej-lMi^i)-  The  emis- 
sion intensity,  x,  and  attenuation  density,  fi,  are  assumed  to  be  deterministic  but 
unknown.  Given  the  emission  data,  y\,y2, . . .  ,yi,  blank  scan  data,  61?  62,  •  •  • ,  bIt 
and  transmission  scan  data,  rur2, . . .  ,rI}  the  goal  is  to  estimate  the  emission 
intensity  and  attenuation  density. 
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It  is  well  known  that  in  many  estimation  problems  ML  estimators  are  consis- 
tent and  asymptotically  unbiased  [44].  Thus,  we  propose  to  estimate  the  unknown 
parameters  using  a  ML  method.  In  other  words,  we  estimate  the  emission  intensity 
and  attenuation  density  by  maximizing  the  log  likelihood  function,  f(x,[i): 

[x*   /»*]  =  arg  max  f(x,fi)  , 
x,  ft  >  0 

where  /(sc,/i)  =  logP[Yi  =  yuY2  =  y2,  ■  ■  ■ ,  Yi  =  y/|as,/i],  and  x*  and  n*  are  ML 

estimates  of  x  and  fi,  respectively.    The  notation  x  >  0  means  that  all  of  the 

elements  of  x  are  nonnegative.  Under  the  model  assumptions,  the  log  likelihood 

function  is  given  by 

f(x,ii)    =    logH  P[Yi  =  yi\x,n]  (3.1) 


t=i 


/     e-l(Px)te-(Lfl),][{px)         {Lfl)t]yt 

=    losII "1 " (3-2) 

i=l  iff 

=    E[-(P*)ie-lLM<  +  M  log(Px)i  -  ViiLrii  -  log(y,-!)]  .   (3.3) 
j=i 

Eq.  (3.1)  follows  from  the  assumption  that  the  data  is  independent.  Note,  since 
the  last  term  in  eq.  (3.3)  does  not  contain  x  or  n,  it  can  be  ignored  in  the 
maximization  procedure.  Thus,  the  above  optimization  problem  is  equivalent  to 
the  following  one: 

(OP)         [x*   (*"]  =  arg  niax  <j){x,n)  , 
where 

<K*iP)  =  Y\-(PxW(L^'  +  yilog(Px)i  -  yiiLn);}  .  (3.4) 

«=1 

The  function  <j){x,fx)  will  be  referred  to  as  the  log  likelihood  function  for  the 

remainder  of  this  dissertation. 
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The  problem  (OP)  is  a  difficult  constrained  optimization  problem  that  must 
be  solved  in  order  to  obtain  the  ML  estimate  of  the  attenuation  density  and  the 
emission  intensity.  In  the  next  section,  we  discuss  some  properties  of  the  log 
likelihood  function  before  presenting  our  method  for  solving  (OP). 


3.2     Concavity  of  The  Log  Likelihood  Function 

Unfortunately,  the  log  likelihood  function,  <f>(x,fi),  is  not  concave  over  the 
feasible  region,  Sx,fi  =  {(x,fi)\x  >  0,  ft  >  0}.  However,  for  a  fixed  nonzero 
vector  z0,  (f>(x,z0)  and  <^>(z0,/i)  are  concave  over  the  sets  Sx  =  {x\x  >  0}  and 
Sfj,  =  {/■*!/*  >  0},  respectively.  These  results  are  proved  below. 

Result  2:  Let  z0  be  a  nonzero  J  x  1  vector.  Then,  <f>(x,  z0)  is  concave 
over  Sx- 

Proof: 

Our  proof  is  based  on  the  well  known  fact  that  a  function  is  concave 
if  its  Hessian  is  negative  semidefinite  [6].  Consider  4>(x,z0),  where 
z0  >  0  is  a  fixed  J  x  1  vector.  It  follows  from  eq.  (3.4)  that  <f>(x,z0) 
is  given  by 

cf>(x,  Zo)  =  J2[-(P*)te-{LZo)t  +  Vi  \og{Px){  -  yi(Lz0)i]  , 
and  its  first  order  partial  derivative  with  respect  to  Xj  equals 


*%*-&&-**"'*■  ™ 
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The  Hessian  of  <f>(x,  z0),  which  we  denote  by  H(x),  is  a  J  x  J  matrix 

whose  (j,  k)th  element  is  given  by 

„,    ,         a      d2<j){x,z0) 

H{x)lk     ~        dx3dxk  ^ 

=     E-^J-   j,*  =  l,2,...,J.  (3.7) 

Consider  the  scalar  mTH{x)m,  where  m  and  x  are  arbitrary  nonzero 
vectors: 

mTH(x)m    =    -EEE^r  (3-8) 

-  S"w-  (3-9) 

Since  for  all  i,  y<  >  0,  (Pm)f  >  0,  and  (Pas),-  >  0,  it  follows  that 
mT H(x)m  <  0,  which  means  that  the  Hessian  is  negative  semidefi- 
nite.  Thus,  <f>(x,  z0)  is  concave  over  Sx- 

Result  3:  Let  z0  be  a  nonzero  J  x  1  vector.  Then,  (f>(z0,  fi)  is  concave 
over  the  set  Su,- 

Proof: 

Similar  to  the  proof  for  Result  2,  we  need  to  prove  that  the  Hessian  of 
<j>(z0,fi)  is  negative  semidefinite  for  z0,  a  nonzero  fixed  J  x  1  vector. 
From  eq.  (3.4),  it  follows  that  <f>(z0,fi)  is  given  by 

<K*o,t*)  =  YH-(Pzo)te-iLW<  +  yilog(Px0)i  -  Vi{L»)i]  , 

t=i 

and  its  first  order  partial  derivative  with  respect  to  fij  equals 

d  l  t  * 

— 4>{z0,n)  =  ^(Pz0)te-(L^L13  -J^ViLij  .  (3.10) 
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The  Hessian  of  <j)(z0,  /a),  which  we  denote  by  H(fi),  is  a  J  x  J  matrix 
whose  (j,  k)th  element  is  given  by 

HMn    =    ^%^  (3.11) 

=    -J2(Pzo),LtJLtke-lLM<,   j,fc  =  l,2,...,J.(3.12) 
«=i 

Consider  the  scalar  mTH(p,)rn,  where  m  and  fx  are  arbitrary  nonzero 

vectors: 

J     J     I 
mTH(ix)m    =     -J£'E'Z(Pzo)iL,jmJLlkmke^L^'  (3.13) 

k=\  j=\  i-\ 

=    -J2(Pzo),(Lm)*e-(LM<.  (3.14) 

!  =  1 

For  all  i,  (Pz0)i  >  0,  e~&t*)i  >  0,  and  {Lm)1  >  0,  hence  it  fol- 
lows that  mTH(n)m  <  0,  which  means  that  the  Hessian  is  negative 
semidefinite.  Thus,  (f>(z0,[i)  is  concave  over  Sfj,. 

3.2.1     Example  to  Illustrate  Concavity  of  the  Log  Likelihood  Function 

To  get  a  sense  of  the  "shape"  of  the  log  likelihood  function,  we  consider  the 
single  voxel  case.  For  this  case,  the  log  likelihood  function  becomes 

<f>(x,  fi)  =  -pxe~lfi  +  y  log(px)  -  yip,   .  (3.15) 

It  is  straightforward  to  show  that  the  Hessian  of  </>(x,p),  which  we  denote  by  if, 
is  given  by 

-3-     Wc-'M 

pie-'"     -piPe-'" 


H 


38 


Consider  the  scalar  mT Hm,  where  m  is  an  arbitrary  nonzero  vector: 

mTHm  =  -\m\  -  pxl2e~lflml  +  2ple-lflmlm2  .  (3.16) 


x2 


From  eq.  (3.16),  it  can  be  seen  that  there  is  no  way  to  guarantee  that  mTHm  is 
less  than  or  equal  to  zero.  For  p  =  1,  /  =  1,  a;  =  40,  p  =  0.049,  y  =  2,  mi  =  0.2167 
and  m2  =  0.01,the  log  likelihood  function  is  not  concave  (i.e.,  mTHm  >  0). 


3.3     Optimization  Procedure 

We  exploit  the  concavity  results  discussed  in  last  section  and  propose  the 
following  optimization  procedure  to  obtain  the  ML  estimates  (i.e.,  solve  (OP)): 

(Step  1)  Given  /xfc_1  and  xk~l ,  determine  xk  by  maximizing  <j)(x,nk~l)  over  5a;- 
In  other  words, 

(PI)        xk  =  argj^  ^(x,^-1). 

(Step  2)  Given  /ifc_1  and  xk,  determine  nk  by  maximizing  <j)(xk,  fi)  over  Sp.  That 
is, 

(P2)         nk  =  m^ct>(xk,n). 

(Step  3)  Iterate  between  Step  1  and  Step  2. 

Step  1  and  Step  2  are  referred  to  as  the  x-step  and  the  p-step,  respectively. 
The  vectors  xk  and  fik  are  the  kth  estimate  of  x  and  fx,  respectively,  where  an 
iteration  consists  of  both  steps  1  and  2.  Unfortunately,  closed  form  solutions  for 
xk  and  fik  are  not  available,  which  means  they  must  be  determined  using  iterative 
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algorithms.  In  the  next  two  sections,  we  develop  algorithms  for  obtaining  xk  and 

..k 


3.3.1     j-step 

Given  fxk~l  >  0,  we  know  that  ^(a;,^-1)  is  concave  over  Sx  from  Result  2. 
Under  the  conditions  in  problem  (PI),  a:  is  a  solution  to  (Pi)  if  and  only  if  the 
Kuhn- Tucker  conditions  [54]  are  satisfied: 


dxj 

d<t>{x,nk-') 


x=x  =  °    if  a?i  ^  0  ,  j  =  1,2, . . . ,  J  , 

\x=x^Q    if  z,  =  0  ,  j  =  1,2,...,  J. 


dxj 
We  can  rewrite  the  above  two  Kuhn-Tucker  conditions  as 

(CI)        Ef^r  =  E^Lflk'i)'^    if*i*0,i-l,2,...IJI 

(C2)        E7^T<E^(jL^"1)'^      if  ^  =  0,i  =  l,2,...,J. 

The  first  condition  (Cl)  suggests  the  following  fixed  point  iteration  for  solving 

(PI): 

(a,)rl  =  (,fc)rJ-MTO.        i  =  1,2,...,J,  (3.17) 


where  (a;*)?1  is  the  mth  estimate  of  the  jth  voxel  of  xk  and  (a!fc)°  =  {xk~l)3  for  all 
ji.  Another  possible  iteration  is 

(**)?+1  =  («*)*^  --££-      Slf    ,':     1,2,..., J. 


^!=1  (P(a;'t)'"), 


However,  in  our  simulations,  this  iteration  "blew  up"  (i.e.,  the  iterates  became 
unbounded).  Consequently,  we  use  the  iteration  in  eq.  (3.17)  to  estimate  x  . 
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The  logic  behind  the  fixed  point  iteration  in  eq.(3.17)  is  discussed  below.  Sup- 
pose (a;fc)™  converges  to  {xk)°°  /  0  ,  then  from  eq.(3.17)  we  have 

(**)r  =  («*)r  w  v-'p^  •  (3-19) 

which  implies  that  condition  (CI)  is  satisfied.  In  other  words,  we  know  that 
the  fixed  point  iteration  in  eq.(3.17)  will  provide  a  solution  satisfying  (CI)  if  it 
converges  to  a  positive  vector.  Additional  motivation  for  the  technique  used  to 
derive  the  algorithm  in  eq.  (3.17)  is  that  the  EMML  algorithm  [28,  46]  and  the 
ISRA  algorithm  [15]  can  be  derived  using  it.  Applying  the  alternating  projection 
method  [10]  to  solve  (Pi)  leads  to  the  same  iterative  equation  as  eq.  (3.17).  This 
iteration  is  guaranteed  to  converge  by  convergence  results  in  Byrne's  paper  [10]. 


3.3.2     //-step 

From  Result  3,  we  know  that  (j)(xh,n)  is  concave  over  Sfi  for  xk  >  0.  Thus, 
/x  is  a  solution  to  (P2)  if  and  only  if  the  following  Kuhn-Tucker  conditions  hold: 

(C3)     0P* V~(XA)i^i  =  I>^i  if^^O,  j  =  l,2,...,J, 

t=i  (=i 

(C4)     £(P*V~(LA)^i  <  I>^i  if/*,=0,  j  =  l,2,...,J. 

t=i  t'=i 

Using  the  same  approach  that  leads  to  the  iteration  in  eq.  (3.17),  we  obtain  the 
following  fixed  point  iteration  for  solving  problem  (P2): 

(/**)?"  =  (M*)rL,'=1      w /        h  ">  i  =  l,2,...,J,  (3.20) 
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where  (ftk)f  is  the  mth  estimate  of  the  ,;th  voxel  of  fik  and  (A*fc)°  =  (/•**  1)j  for  all 
j.  Another  possible  iterative  equation  is 

(M>r,  =  (Ay__^^__,  J  =  l,2,...,J.  (3.21) 

In  our  simulations,  this  iteration  "blew  up".  Consequently,  we  estimate  fik  using 
the  iteration  in  eq.  (3.20). 

In  the  next  section,  we  investigate  the  behavior  of  the  log  likelihood  function 
as  a  function  of  the  iterates. 


3.3.3     Monotonically  Increasing  Value  of  the  Log  Likelihood  Function  in  j-step 

From  the  log  likelihood  function  given  in  eq.(3.4)  and  the  iterative  equation  in 
eq.(3.17),  we  can  prove  that  the  log  likelihood  function  increases  monotonically 
with  successive  iteration  in  the  a>step.  Unfortunately,  this  result  cannot  be  proven 
for  the  /x-step. 

Result  5:  Let  z0  be  a  fixed  nonzero  J  x  1  vector.  Then, 

<f>((xk)m+\zo)-<l>((xkr,zo)>0. 
Proof: 
From  eq.(3.4),  it  is  easy  to  see  that 

ct>{{xkr^\z0)-<j>{{xkr,z0) 

=  U-(P(*T+1)i  +  (P(«*rM«i  +  E n log \p{xl)m))t' . 

where  a,-  =  e~(       °''. 
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Hence,  using  the  iterative  equation  in  eq.(3.17),  it  follows  that  the 
second  term  of  the  above  derivation  can  be  written  as 

|>'°«  wn   -  £Klogg(P(^)7 

/  J         p 


=  Zy^nj^yr^ 


I  J     P-(xk)m 


TL, 


ViPi 


where  u^  =  — _/     ^ — *— "-. 

p.    /jg/e\m 

Let  /?t-j  =    p     ggr  ,  £i=i  Aj  —  1-  Using  the  fact  that  the  log  function 
is  concave,  we  have 

E  Vi  log  E  fa")  >  E  V>  E  Ai  log  «?  •  (3-22) 

1=1  j=l  1=1     j=l 

Multiplying  the  right  side  of  eq.(3.22)  by  ,^'r' "'  '  .,  the  value  of  right 

t\^I  p.    \ 

side  of  eq.(3.22)  does  not  change  since  fo=1  °'  'J,  =  1-  Then,  it  follows 
that 


Ey.E^iog^  =  Zyi(£<*iPii)Zh^T' 

i=l  j  =  l  !  =  1  !=1  l^i  =  l  at*>j 

j=i     t=i  ^,=i  ««^tj 

=  E4(E^-K*iog^ 

j=l  i=l 
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j-i    t=i 


=  E^E(j&:-E<*E*i<«T 


i=l  t'=l 

Since  E]L,  2/*  log  £/=1  ftf"?  >  £/=1  Sfc  E/=,  A,  logt^,  we  have 

7  (P(rk}m+1\-         '  l 

X>'°6    (p"J)m),    £  g>  -  E".(P(^D.  ■  (3-23) 

Then,  we  get 

>  EL  w  -  ELKfPf^D.  +  «,(P(*T)<  -  at(P(xkr+i)l] 

From  eq.(3.4),  we  know  that 

gaflVr**). = E  "l^l')"'  ■  g»  ■        <3-24> 

From  eq.(3.24),  we  obtain 

E*-i:o»(P(«*)-+1)i«0  (3.25) 

,=i  i=i 

Thus,  we  get  the  desired  result. 

^((a5fer+i,zo)-^((a!fcr,zo)>o. 

In  the  next  section,  an  interesting  property  of  the  ML  estimator  is  presented. 
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3.3.4     Property  of  the  ML  Estimator 

Result  4-'  Suppose  there  exists  a  solution  (x*,fi*)  to  (OP).  Then, 

ViPi, 


ZUPx*)te-^^LtJ 


<  1  ,  (3.26) 

<  1,   7  =  1,2,..., J.    (3.27) 


E,=i  ViLij 
Proof: 
Consider  the  log  likelihood  function 

tf«,/i)  =  E[-(P*)iC-^^  +yilog(P*)i-yi(Lp)i] 


!  =  1 


Suppose  that  there  exists  an  optimal  solution  (x*,fi*)  and  let  t  be 
a  nonnegative  scalar,  and  x  and  fi  be  nonnegative  vectors.  Then, 
<j)(x*  +  tx,  n*  +  tfi)  <  <j)(x*,  fi*).  From  the  above  relationship,  we  get 

-[(Px')i  +  *(Pz).]e~t(jL/X*)'+<(Z'At)']  +  (Px*)^L^}  <  0 
Since  (Px*)^1^'  >  0  for  all  t,  we  have 

zLivd<*[{P*;)%$x)i] 

-  ty,(Ln)i  -  [(Px*){  +  t(Px)i]e^L^'+^L^}  <  0         (3.28) 
Dividing  eq.(3.28)  by  t,  we  get 

i:{y^og[i+^^yy^^  <  o 

(3.29) 
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Now  letting  t  — >  0  and  using  the  L'Hoptial's  rule,  we  obtain 

Eij^V  -  ViiWi  +  (P*)KLriie-<LW  -  (Px),e^L^}  <  0. 

(3.30) 
Since  eq.(3.30)  must  hold  for  any  \i  >  0  and  x  >  0,  we  consider  the 
case  where  x  —  ej  (ej  is  a  J  x  1  vector  of  zeros  except  for  the  jth 
element  which  is  a  one)  and  fj,  =  0.  For  this  case,  eq.  (3.30)  becomes 

h^V-^Lfit)'^]<^   J  =  1,2,..., J.  (3.31) 

i=l    \^X    )i 


Similarly,  for  x  =  0  and  \x  =  ej,  eq.  (3.30)  yields 

Y,[-ViLij  +  (Px*)^^'^}  <  0,   j  =  1, 2, . . . ,  J  .        (3.32) 


i=i 


Note,  eq.  (3.31)  and  (3.32)  are  only  necessary  conditions  for  (x*,[i*) 
to  be  a  maximum  of  (j)(x,  fi)  over  the  feasible  region.  Since  x*,  fi*  >  0, 
it  follows  from  eq.  (3.31)  and  eq.  (3.32)  that 


Tl        y,Pt] 

i-Jl~\     (  Pj>.1 

<1,   j  =  l,2,...,J,  (3.33) 


^e-iLfhPij 
Ul7fJ'     , ^<l,   ;  =  1,2,...,J.  (3.34) 

Ei=i  y.^«i 

In  the  next  chapter,  we  introduce  a  new  least  squares  reprojection  attenuation 
correction  method. 


CHAPTER  4 
LEAST  SQUARES  REPROJECTION  METHOD 

In  this  chapter,  we  introduce  a  reprojection  method,  where  the  attenuation 
density  is  estimated  using  a  least  squares  (LS)  algorithm.  The  algorithm  is  an  ex- 
tension of  one  proposed  by  Anderson,  et.  al.  [4]  for  X-ray  computed  tomography. 

4.1     LS  Algorithm  for  Estimating  the  Attenuation  Density 

As  discussed  in  Section  2.4  and  2.5,  let  6t-  denote  the  number  of  photons  de- 
tected in  tube  i  during  the  blank  scan.  It  is  assumed  that  6,  is  an  observation  of 
a  random  variable,  1?,-,  where  B{  is  Poisson  with  mean  [/,  >  0.  It  is  also  assumed 
that  the  number  of  photons  detected  in  tube  i  during  the  transmission  scan,  rt-,  is 
an  observation  of  a  random  variable,  i?,,  where  i?,  is  Poisson  with  mean  Uia^111) 
(recall  a,  =  e~(  M').  Given  the  observed  blank  scan  data,  b\,  b2, . . . ,  bj,  and  the 
observed  transmission  data,  ri,r2, . . . , r/,  the  goal  is  to  estimate  the  attenuation 
density,  \i. 

It  is  straightforward  to  show  that  6;  and  r,-  are  the  ML  estimates  of  £/,•  and 
UiCti^),  respectively.  Consequently,  6,  «  £/,-  and  r,-  «  C/.-a.-^).  These  approxi- 
mations lead  to  the  following  one:  a,-  «  ^i2^)-  Taking  the  logarithm  on  both  sides 
of  this  relationship  and  multiplying  the  result  by  negative  one,  we  get  (Lfi)t  «  d,, 
where  d{  =  —  log(^i(^!l)).  Exploiting  this  observation,  we  propose  to  estimate  the 
attenuation  density  by  minimizing  the  2-norm  distance  between  d  and  (Lfi): 


(P3)         //*  ^mii^  p(p)  , 


46 


47 


where  p(fi)  =  E.L^d,-  — (I/fx),-)2,  A**  iS  the  LS  estimate  of  fj,,  and  d  =  [di,d2,.  ..,di]. 
It  is  straightforward  to  show  that  p(fi)  is  convex  over  Sp. 

A  vector  ft  is  a  solution  to  (P3)  if  and  only  if  it  satisfies  the  following  Kuhn- 
Tucker  conditions: 

(C5)        "£Ll3(di-(Lii)t)  =  0  if  ^  ^  0  ,  i  =  1,2, ...,  J  , 

(ce)      £-Lo-(4-(L£).-)>o  if  ^  =  o,j  =  i,2,...,y. 

i=i 

As  before  (see  Chapter  3),  the  first  condition  (C5)  motivates  the  following  fixed 
point  iteration: 

where  //^  is  the  kth  estimate  of  fij.  This  iteration  has  exactly  the  same  form  as 
the  iterative  image  space  reconstruction  algorithm  for  PET  [15].  Since  it  has  been 
proven  that  the  iterative  image  space  reconstruction  algorithm  converges  [40],  the 
proposed  LS  algorithm  converges  as  well.  In  the  next  section,  we  discuss  how  the 
LS  estimate  of  the  attenuation  density  is  used  to  correct  for  attenuation. 

4.2     LS  Reprojection  Method 

Given  the  blank  scan  data  and  the  transmission  scan  data,  we  estimate  the 
attenuation  density  using  eq.  (4.1).  Let  ft  be  the  resulting  estimate.  Then,  the 
survival  probabilities  are  estimated  by  reprojecting  ft: 

a,  =  e-£/=iL'^    ,t=i  1,2,..., I.  (4.2) 

Next,  the  estimated  survival  probabilities  are  used  to  create  a  modified  probability 
matrix,  P,  where  P,j  =  a,'P,j.  The  emission  intensity  is  then  estimated  by  using 
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the  EMML  algorithm  with  the  modified  probability  matrix: 

It  is  straightforward  to  show  that  this  iteration  is  equivalent  to 

*?' =**,''■   ^V^   i=lA...,J.  (4.4) 

In  the  next  chapter,  simulations  are  used  to  compare  the  performance  of  the 
ML  method,  LS  reprojection  method,  and  EMML  algorithm  with  corrected  data. 


CHAPTER  5 
SIMULATION 

Using  synthetic  data  and  real  data,  we  assess  the  proposed  methods  and  com- 
pare them  to  the  EMML  algorithm  with  corrected  data.  Dr.  John  Votaw,  As- 
sociate Professor  of  Radiology  at  the  Emory  University  Hospital  provided  the 
real  data.  This  data  was  collected  on  the  Siemens-CTI  ECAT  EXACT  921  PET 
scanner.  The  probabilities,  P,j,  and  line  length,  L{j,  were  determined  from  the 
diameter  of  the  detector  ring,  number  of  detectors,  the  size  of  ROI  (the  ROI  is 
always  a  square  centered  in  the  middle  of  the  detector  ring),  and  the  number  of 
voxels.  Specially,  the  probabilities  are  calculated  using  the  angle-of-view  method 
[46].  It  should  be  noted  that  the  probabilities  and  line  lengths  are  calculated  un- 
der the  assumption  that  the  detectors  are  contiguous.  The  number  of  detectors 
for  the  simulated  scanner  and  the  Siemens-CTI  ECAT  EXACT  921  scanner  are 
300  and  384,  respectively.  The  computational  expense  of  an  iteration  for  the  ML 
method,  LS  reprojection  method,  and  the  EMML  algorithm  with  corrected  data 
are  essentially  the  same.  Consequently,  the  computational  complexity  of  the  al- 
gorithms can  be  compared  by  considering  the  total  number  of  iterations  required 
by  each  one. 
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5.1     Synthetic  Data 

For  the  synthetic  data,  we  use  a  cardiac  emission  phantom  (see  Fig.  5.3  (a)), 
chest  emission  phantom  (see  Fig.  5.7  (a)),  cardiac  attenuation  phantom  (see  Fig. 
5.5  (a)),  and  chest  attenuation  phantom  (see  Fig.  5.8  (a)).  All  of  the  phantoms  are 
of  size  100  x  100,  and  the  voxel  size  for  the  cardiac  phantoms  and  chest  phantoms 
are  0.422  cm  and  0.437  cm,  respectively.  The  lungs,  spine,  and  soft  tissue  region 
in  the  cardiac  attenuation  phantom  have  densities  0.025  cm-1,  0.16  cm-1,  and 
0.096  cm-1,  respectively.  The  low  density  regions  (light  color  region)  and  the 
high  density  regions  (dark  color  region)  in  the  chest  attenuation  phantom  have 
densities  0.035  cm  and  0.095cm_1,  respectively.  The  diameter  of  the  detector 
ring  for  the  cardiac  phantoms  and  chest  phantoms  is  59.68  cm  and  61.8  cm, 
respectively.  The  ROI  for  the  cardiac  phantoms  and  chest  phantoms  is  42.2  cm 
and  43.7  cm,  respectively.  The  number  of  tubes  is  /  =  44850. 

We  generate  the  emission  scan,  transmission  scan,  and  blank  scan  data  in  the 
following  manner.  The  observed  coincidence  for  the  ith  tube,  ?/,,  is  the  observation 
of  a  Poisson  random  variable  with  mean  (Pa^-a,  (recall  a,  =  e~(  M*).  The 
emission  phantoms  and  attenuation  phantoms  specify  x  and  ft,  respectively.  The 
observed  count  for  the  blank  scan  for  the  iih  tube,  6,-,  and  the  transmission  scan 
for  the  ith  tube,  r,,  are  Poisson  random  variables  with  means  £/,-  and  £/,•<*; (^J, 
respectively.  We  specify  the  total  count  of  the  transmission  scan  to  be  R.  and  use 


^  "  (H-i ")(*)/' *-1'2"-"7" 
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5.1.1  ML  Method  Implementation  of  Algorithms 

The  term  local  iteration  refers  to  the  iterations  used  to  determine  xk  and  fik 
in  Step  1  (x-step)  and  Step  2  (//-step),  respectively  (see  (PI)  and  (P2)  in  Section 
3.3).  The  term  global  iteration  refers  to  a  single  execution  of  both  the  x-step  and 
the  //-step.  For  the  simulations  using  synthetic  data,  we  run  8  local  iterations  for 
the  x-step,  5  local  iterations  for  the  /f-step,  and  2  global  iterations.  Thus,  a  total 
of  21  iterations  and  26  iterations  are  used  to  obtain  emission  intensity  (i.e.,  x2) 
and  attenuation  density  (i.e.,  ft2),  respectively. 

The  ML  method  is  initialized  using  (aj0)^  =  £i=i  ft,  and  (/*°)i  =  0-0214> 
except  that  e'^V-0)*  is  replaced  with  a°  =  ^(^).  As  state  in  Section  3.3.1  and 
Section  3.3.2,  (**)?  =  {xk-x)j  and  (ji*)J  -  (/i*"1);- 

5.1.2  EMML  Algorithm  with  Corrected  Data 

For  the  EMML  algorithm  with  corrected  data,  the  iterative  equation  is  given 
below: 


k  (p*k)i ' 


•?■- -•}£?£*:.  ^ 


where  y  =  fS    It  is  initialized  as  follows:    x?  =   ^'t1  y'  for  all  j.     We  run  20 

*  ati  3  1* 

iterations  of  the  EMML  algorithm  with  corrected  data. 

5.1.3     LS  Reprojection  Method 

To  initialize  the  iteration  in  eq.  (4.1),  we  use  a  uniform  initial  estimate:  p,°  = 
0.0214  for  all  j.  We  run  20  iterations  to  get  the  LS  estimate  of  ft,  ft.  Then,  we 
run  20  iterations  of  the  algorithm  shown  in  eq.  (4.4).  Thus,  the  LS  reprojection 
method  requires  40  iterations  to  get  an  emission  image. 
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The  total  emission  count,  typical  transmission  scan  count,  and  short  transmis- 
sion scan  count  are  1.5  million,  0.5  million  and  0.125  million,  respectively.  The 
blank  scan,duration,  the  typical  transmission  scan  duration,  and  the  short  trans- 
mission scan  duration  are  60  minutes,  20  minutes,  and  5  minutes,  respectively.  We 
simulate  the  ML  method,  LS  reprojection  method,  and  the  EMML  algorithm  with 
corrected  data  using  typical  transmission  scan  data  with  and  without  smoothing 
and  the  short  transmission  scan  data  with  and  without  smoothing.  In  the  next 
section,  we  show  the  results  using  non-smoothed  typical  transmission  scan  data. 


5.1.4     Non-smoothed  Typical  Transmission  Scan  Data 

We  present  the  results  using  a  non-smoothed  typical  transmission  scan  data  in 
this  section.  The  emission  image  of  the  cardiac  emission  phantom  using  the  ML 
method  after  the  1st  global  iteration,  2nd  global  iteration,  and  3rd  global  iteration 
are  shown  in  Fig.  5.1  (b),  Fig.  5.1  (c),  and  Fig.  5.1  (d),  respectively.  It  can  be 
seen  that  the  emission  images  improve  with  increasing  iteration.  However,  after 
3  global  iterations,  the  emission  images  become  noisy  and  less  accurate,  which 
indicates  that  a  stopping  rule  needs  to  be  found.  This  problem  is  a  common  one 
for  statistically  based  algorithms. 

The  reconstructions  of  the  cardiac  attenuation  phantom  using  the  ML  method 
after  the  1st  global  iteration,  2nd  global  iteration,  and  3rd  global  iteration  are  shown 
in  Fig.  5.2  (b),  Fig.  5.2  (c),  and  Fig.  5.2  (d),  respectively.  In  the  attenuation 
images,  the  lungs  are  not  visible,  and  only  the  contour  and  the  spine  can  be  seen. 
Also,  there  are  artifacts  outside  the  ROI. 
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Figure  5.1.  Reconstruction  of  cardiac  emission  phantom  using  non-smoothed  typ- 
ical transmission  scan  data  (a)  phantom  (b)  1st  global  iteration  (c)  2nd  global 
iteration  (d)  3rd  global  iteration 
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Figure  5.2.  Reconstruction  of  cardiac  attenuation  phantom  using  non-smoothed 
typical  transmission  scan  data  (a)  attenuation  phantom  (b)  1st  global  iteration 
(c)  2nd  global  iteration  (d)  3rd  global  iteration 
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For  the  cardiac  emission  phantom,  the  results  of  the  EMML  algorithm  with 
corrected  data,  the  LS  reprojection  method,  and  the  ML  method  are  shown  in 
Fig.  5.3  (a),  Fig.  5.3  (b),  Fig.  5.3  (c),  and  Fig.  5.3(d),  respectively.  From 
this  figure,  it  can  be  seen  that  the  ML  method  yields  an  image  with  better  con- 
trast and  resolution  than  the  EMML  algorithm  with  corrected  data  and  the  LS 
reprojection  method  (note  the  lungs  and  the  heart  at  the  top  center  of  the  cardiac 
phantom).  It  also  can  be  seen  that  the  EMML  algorithm  with  corrected  data 
does  not  successfully  resolve  either  the  lungs  or  the  heart.  In  the  emission  image 
of  the  LS  reprojection  method,  the  spine  can  be  seen  clearly.  However,  the  LS 
reprojection  method  also  provides  incorrect  information  as  seen  by  the  fact  that 
the  lungs  have  a  higher  activity  concentration  than  the  soft  tissues  around  them 
(see  the  phantom  in  Fig.  5.3  (a)). 

In  Fig.  5.4,  we  show  the  line  plots  of  the  50th  row  of  the  emission  images  in  Fig. 
5.3  .  We  observe  that  both  the  LS  reprojection  method  and  the  EMML  algorithm 
with  corrected  data  are  more  biased  than  the  ML  method  (note  the  plot  between 
52  and  62,  which  is  inside  the  right  lung  area).  From  this  line  plot,  we  find  that 
the  ML  method  not  only  produces  emission  images  with  better  contrast  but  also 
gives  a  more  accurate  estimate  of  the  emission  intensity  than  other  methods. 

The  EMML  algorithm  with  corrected  data  does  not  reconstruct  attenuation 
images,  so  we  only  show  the  attenuation  images  of  the  ML  method  and  the  LS 
reprojection  method  in  Fig.  5.5  (b)  and  Fig.  5.5  (c),  respectively.  The  ideal 
attenuation  phantom  is  shown  in  Fig.  5.5  (a).  We  can  see  that  the  attenuation 
density  reconstruction  of  the  ML  method  in  Fig.  5.5  (b)  has  a  lot  of  artifacts  and 


56 


100 


20  40  60  80  100 

(a) 


100 


100 


Figure  5.3.  Reconstruction  of  cardiac  emission  phantom  using  non-smoothed  typ- 
ical transmission  scan  data  (a)  emission  phantom  (b)  EMML  algorithm  with  cor- 
rected data  (c)  LS  reprojection  method  (d)  ML  method 
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Dot:  EM-Standard  Correction  Method 
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Figure  5.4.    Line  plot  of  the  50th  row  of  the  reconstruction  of  cardiac  emission 
phantom  using  non-smoothed  typical  transmission  scan  data 
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does  not  resolve  the  lungs,  spine,  and  soft  tissue.  The  LS  reprojection  method 
resolves  lungs  and  soft  tissues  but  also  has  artifacts. 

We  also  show  the  emission  images  of  the  ML  method  initialized  with  a  LS 
estimate  of  the  attenuation  density.  Since  the  fi-step  is  supposed  to  determine  a 
better  estimate  of  /x,  we  want  to  check  if  a  better  initial  fi  will  correct  the  artifact 
problem  of  the  attenuation  images  and  provide  better  emission  images  for  the  ML 
method.  In  Fig.  5.6  (b),  we  can  observe  that  spine  is  clearer,  but  the  lung  area  is 
incorrect.  However,  the  attenuation  image  in  Fig.  5.6  (d)  is  better  than  Fig.  5.5 
(b)  (resolve  lungs,  soft  tissue  but  no  spine). 

The  results  of  the  EMML  algorithm  with  corrected  data,  the  LS  reprojection 
method  and  the  ML  method  using  the  chest  phantom  are  shown  in  Fig.  5.7  (b), 
5.7  (c),  and  5.7  (d),  respectively.  From  this  figure,  we  can  observe  that  the  ML 
method  and  the  LS  reprojection  method  produce  comparable  emission  images. 
Both  images  have  better  contrast  and  resolution  than  the  EMML  algorithm  with 
corrected  data.  Note,  the  emission  image  produced  by  the  EM  algorithm  with 
corrected  data  is  very  noisy. 

We  present  the  attenuation  images  of  the  ML  method  and  the  LS  reprojection 
method  for  the  chest  phantom  in  Fig.  5.8  (b)  and  Fig.  5.8  (c),  respectively.  The 
ideal  attenuation  phantom  is  shown  in  Fig.  5.8  (a).  For  the  chest  phantom,  the 
ML  method  provides  better  contrast  and  more  details  (big  and  small  spots  in  the 
center)  than  the  LS  reprojection  method.  The  LS  reprojection  method  fails  to 
resolve  any  information. 
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Figure  5.5.  Reconstruction  of  cardiac  attenuation  phantom  using  non-smoothed 
typical  transmission  scan  data  (a)  attenuation  phantom  (b)  ML  method  (c)  LS 
reprojection  method 
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Figure  5.6.  ML  method  (a)  emission  phantom  (b)  emission  image  (c)  attenuation 
phantom  (d)  attenuation  image 
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Figure  5.7.  Reconstruction  of  chest  emission  phantom  using  non-smoothed  typical 
transmission  scan  data  (a)  emission  phantom  (b)  EMML  algorithm  with  corrected 
data  (c)  LS  reprojection  method  (d)  ML  method 
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Figure  5.8.  Reconstruction  of  chest  attenuation  phantom  using  non-smoothed 
typical  transmission  scan  data  (a)  attenuation  phantom  (b)  ML  method  (c)  LS 
reprojection  method 
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In  practice,  the  transmission  scan  data  is  smoothed  even  though  it  is  ad  hoc  to 
do  so  [34].  In  next  section,  we  show  the  results  of  the  simulations  using  smoothed 
transmission  scan  data. 


5.1.5     Smoothed  Typical  Transmission  Scan  Data 

Since  the  survival  probability  estimated  by  the  ratio  of  transmission  scan  data 
and  blank  scan  data  is  usually  noisy  and  biased,  a  popular  way  to  address  this 
problem  is  to  smooth  the  transmission  scan  data  [13,  34,  41].  We  use  a  1  x  3  box 
car  filter,  [|  |  |]  to  smooth  the  transmission  scan  data. 

We  show  the  emission  images  of  the  EMML  algorithm  with  corrected  data, 
the  LS  reprojection  method  and  the  ML  method  for  the  cardiac  phantom  using 
smoothed  typical  transmission  scan  data  in  Fig.  5.9  (b),  5.9  (c),  and  Fig.  5.9  (d), 
respectively.  Note  that  the  LS  reprojection  method  resolves  more  details  of  the 
lungs,  heart  and  the  spine  than  the  ML  method.  Both  methods  produce  images 
much  better  than  the  EMML  algorithm  with  corrected  data. 

This  study  demonstrates  that  smoothing  the  transmission  scan  data  improves 
the  quality  of  the  emission  images.  In  order  to  compare  the  accuracy  of  the 
emission  image,  in  Fig.  5.10,  we  show  the  line  plots  of  the  50th  row  of  the  images 
in  Fig.  5.9.  We  can  see  that  the  line  plot  of  the  LS  reprojection  method  and 
the  EMML  algorithm  with  corrected  data  both  provide  worse  estimates  of  the 
emission  intensity  than  the  ML  method. 

We  also  compare  the  attenuation  images.  We  show  the  result  of  the  ML  method 
and  the  LS  reprojection  method  in  Fig.  5.11  (b)  and  Fig.  5.11  (c),  respectively. 
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Figure  5.9.  Reconstruction  of  cardiac  emission  phantom  using  smoothed  typical 
transmission  scan  data  (a)  emission  phantom  (b)  EMML  algorithm  with  corrected 
data  (c)  LS  reprojection  method  (d)  ML  method 
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Figure  5.10.  Line  plot  of  the  50"1  row  of  reconstruction  of  cardiac  emission  phan- 
tom using  the  smoothed  transmission  scan  data 
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The  ML  method  still  produces  an  attenuation  image  with  heavy  artifacts  while 
the  LS  method  provides  an  image  with  better  contrast  and  quality.  Compared  to 
the  attenuation  images  using  non-smoothed  transmission  scan  data  in  Fig.  5.5, 
both  attenuation  images  show  more  features. 

The  emission  and  attenuation  image  of  the  ML  method  using  an  initial  esti- 
mate for  fi  was  obtained  from  the  LS  reprojection  method  are  shown  in  Fig.  5.12 
(b),  Fig.  5.12  (d),  respectively.  The  ML  method  produces  emission  image  with 
better  contrast  in  the  heart  area  in  than  Fig.  5.9  (b).  The  attenuation  image 
provides  more  information  than  both  attenuation  images  in  Fig.  5.11  (note  the 
black  spine  at  the  center  below). 

The  emission  images  of  the  EMML  algorithm  with  corrected  data,  the  LS 
reprojection  method  and  the  ML  method  of  the  chest  phantom  using  smoothed 
typical  transmission  scan  data  are  shown  in  Fig.  5.13  (b),  5.13  (c),  and  5.13  (d), 
respectively.  We  can  see  that  the  LS  reprojection  method  and  the  ML  method 
produce  comparable  images.  Both  methods  provide  better  contrast  and  the  image 
quality  than  the  EMML  algorithm  with  corrected  data. 

The  attenuation  images  of  the  ML  method  and  the  LS  reprojection  method 
using  the  typical  smoothed  transmission  data  for  chest  attenuation  phantom  are 
shown  in  Fig.  5.14  (b)  and  Fig.  5.14  (c),  respectively.  The  ML  method  resolves 
lungs,  a  heart  and  soft  tissue  while  the  LS  method  fails  to  provide  any  information. 
The  attenuation  image  of  the  ML  method  is  smoother  than  Fig.  5.8  (b)  but  the 
contrast  and  quality  do  not  improve. 
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Figure  5.11.  Reconstruction  of  cardiac  attenuation  phantom  using  the  smoothed 
typical  transmission  scan  data  (a)  attenuation  phantom  (b)  ML  method  (c)  LS 
reprojection  method 
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Figure  5.12.  ML  method  (a)  emission  phantom  (b)  emission  image  (c)  attenuation 
phantom  (d)  attenuation  image 
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Figure  5.13.  Emission  image  of  the  chest  phantom  using  smoothed  typical  trans- 
mission scan  data  (a)  emission  phantom  (b)  EMML  algorithm  with  corrected  data 
(c)  LS  reprojection  method  (d)  ML  method 
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Figure  5.14.  Reconstruction  of  chest  attenuation  phantom  using  smoothed  typical 
transmission  scan  data  (a)  attenuation  phantom  (b)  ML  method  (c)  LS  reprojec- 
tion  method 
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In  order  to  compare  the  methods  numerically,  we  compute  their  peak-signal- 
to-noise-ratio  (PSNR)  [43]  for  the  different  images.  PSNR  is  defined  by 

PSNR(dB)  =  101oglo^,  (5.2) 

mse 

where  xpeak  is  the  peak  value  of  the  phantom,  mse  —  jY1j=i{xj  -  £j)2>  xj  1S 

the  jth  voxel  in  phantom,  and  Xj  is  an  estimate  of  Xj.   We  compare  the  PSNR 

of  the  methods  using  the  typical  transmission  scan  data  in  Table  5.1.      We  can 


clB 

cardiac  phantom 

chest  phantom 

Algorithm 

Standard 

ML 

LS-Rep. 

Standard 

ML 

LS-Rep. 

Non-smooth 

-7.5 

12.51 

-5.92 

-6.3 

12.73 

-4.8 

Smooth 

-15.5 

8.62 

-5.62 

-6.7 

12.78 

-4.56 

Table  5.1.  PSNR  (dB)  of  emission  images  using  typical  transmission  scan  data 

observe  that  the  emission  images  of  the  ML  method  have  the  highest  PSNR  and 
the  EMML  algorithm  with  corrected  data  has  the  lowest  PSNR. 

Another  application  of  the  proposed  methods  is  the  case  of  short  duration 
transmission  scan  data.  In  the  next  two  sections,  we  show  the  results  using  a  5 
minutes  short  transmission  scan  data  with  and  without  smoothing. 


5.1.6     Non-smoothed  Short  Transmission  Scan  Data 

In  this  section,  we  show  the  results  using  a  5  minutes  transmission  scan.  The 
emission  images  of  the  EMML  algorithm  with  corrected  data,  the  LS  reprojection 
method  and  the  ML  method  of  cardiac  emission  phantom  using  non-smoothed 
typical  transmission  scan  data  are  shown  in  Fig.  5.15  (b),  5.15  (c),  and  5.15  (d), 


T> 


respectively.  From  this  figure,  we  know  that  the  ML  method  resolves  the  heart 
and  lungs  better  than  the  the  EMML  algorithm  with  corrected  data  and  the  LS 
reprojection  method.  The  emission  image  of  the  LS  reprojection  method  produces 
incorrect  information  in  the  lungs  region. 

We  also  compare  the  attenuation  images.  We  show  the  result  of  the  ML  method 
and  the  LS  reprojection  method  in  Fig.  5.16  (b)  and  Fig.  5.16  (c),  respectively. 
We  can  see  that  the  attenuation  image  of  ML  method  in  Fig.  5.16  (b)  has  a  lot  of 
artifacts  and  only  resolves  spine  (black  spot  at  the  center  below).  The  attenuation 
image  of  the  LS  reprojection  method  in  Fig.  5.16  (c)  shows  the  basic  contours 
better  than  the  ML  method  but  also  suffers  from  the  artifacts.  Both  methods 
provide  worse  attenuation  images  when  using  short  transmission  scan  data. 

The  emission  and  attenuation  image  of  the  ML  method  initialized  using  the 
LS  estimate  of  fi  are  shown  in  Fig.  5.17  (b),  Fig.  5.17  (d),  respectively.  The  ML 
method  resolves  the  heart  area  better  than  Fig.  5.15  (d)  but  produces  the  wrong 
information  of  the  lungs.  The  attenuation  image  provides  more  information  and 
less  artifacts  than  Fig.  5.16  (c). 

The  results  of  the  EMML  algorithm  with  corrected  data,  the  LS  reprojection 
method  and  the  ML  method  of  chest  emission  phantom  using  non-smoothed  short 
transmission  scan  data  are  shown  in  Fig.  5.18  (b),  5.18  (c),  and  5.18  (d),  respec- 
tively. From  this  figure,  we  can  observe  that  the  ML  method  that  produces 
image  with  better  contrast,  resolution  than  the  EMML  algorithm  with  corrected 
data  and  the  LS  reprojection  method.  The  EMML  algorithm  with  corrected  data 
and  LS  reprojection  method  both  produce  image  with  wrong  information.    The 
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Figure  5.15.  Reconstruction  of  cardiac  emission  phantom  using  the  non-smoothed 
short  transmission  scan  data  (a)  emission  phantom  (b)  EMML  algorithm  with 
corrected  data  (c)  LS  reprojection  method  (d)  ML  method 
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Figure  5.16.  Reconstruction  of  cardiac  attenuation  phantom  using  non-smoothed 
short  transmission  scan  data  (a)  attenuation  phantom  (b)  ML  method  (c)  LS 
reprojection  method 
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Figure  5.17.  ML  method  (a)  emission  phantom  (b)  emission  image  (c)  attenuation 
phantom  (d)  attenuation  image 


76 


100 


20 
40 
60 
80 
100 


7-5 


20  40  60  80  100 

(C) 


100 


100 


100 


100 


Figure  5.18.  Reconstruction  of  chest  emission  phantom  using  non-smoothed  short 
transmission  scan  data  (a)  emission  phantom  (b)  EMML  algorithm  with  corrected 
data  (c)  LS  reprojection  method  (d)  ML  method 
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emission  image  of  the  LS  reprojection  method  is  smoother  than  the  one  of  the 
EMML  algorithm  with  corrected  data. 

We  present  the  attenuation  images  of  the  ML  method  and  the  LS  reprojection 
method  in  Fig.  5.19  (b)  and  Fig.  5.19  (c),  respectively.  The  LS  method  does  not 
provide  too  much  information  while  the  ML  method  provides  a  noisy  attenuation 
image.  However,  the  ML  method  resolves  the  high  density  region  at  the  center. 

In  the  next  section,  we  show  the  results  of  the  simulations  using  smoothed 
short  transmission  scan  data. 


5.1.7     Smoothed  Short  Transmission  Scan  Data 

In  this  section,  we  show  the  results  with  short  transmission  scan  smoothed 
by  the  same  1x3  box  car  filter,  [|  |  |].  The  emission  images  of  the  EMML 
algorithm  with  corrected  data,  the  LS  reprojection  method  nd  the  ML  method  of 
cardiac  phantom  using  smoothed  short  transmission  scan  data  are  shown  in  Fig. 
5.20  (b),  5.20  (c),  and  5.20  (d),  respectively.  Again,  we  observe  that  We  can  see 
that  the  LS  reprojection  method  provides  the  comparable  emission  image  with 
the  ML  method  in  the  sense  of  the  contrast  and  image  quality.  Note  that  the 
LS  reprojection  method  resolves  more  details  of  the  lungs,  heart  and  spine  but 
also  suffers  from  quantitative  inaccuracy.  Especially,  the  LS  reprojection  method 
provides  clearer  heart  area  than  other  methods.  The  ML  method  and  the  LS 
reprojection  method  both  methods  produce  images  much  better  than  the  EMML 
algorithm  with  corrected  data. 
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Figure  5.19.  Reconstruction  of  chest  phantom  using  non-smoothed  short  trans- 
mission scan  data  (a)  attenuation  phantom  (b)  ML  method  (c)  LS  reprojection 
method 
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Figure  5.20.  Reconstruction  of  cardiac  emission  phantom  using  smoothed  short 
transmission  scan  data  (a)  emission  phantom  (b)  EMML  algorithm  with  corrected 
data  (c)  LS  reprojection  method  (d)  ML  method 
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We  show  the  attenuation  images  of  the  ML  method  and  the  LS  reprojection 
method  in  Fig.  5.21  (b)  and  Fig.  5.21  (c),  respectively.  The  attenuation  image 
of  the  LS  method  provides  the  shape  of  the  attenuation  map  while  the  image  of 
the  ML  method  suffers  from  heavy  artifacts. 

The  emission  and  attenuation  image  of  the  ML  method  initialized  with  ft 
are  shown  in  Fig.  5.22  (b),  Fig.  5.22  (d),  respectively.  The  ML  method  almost 
produces  the  best  emission  images  for  the  all  short  transmission  scan  experiments. 
In  particularly,  it  provides  a  very  clear  reconstruction  of  the  heart  area  (see  Fig. 
5.9  (b)). 

The  emission  images  of  the  EMML  algorithm  with  corrected  data,  the  LS 
reprojection  method  and  the  ML  method  of  the  chest  emission  phantom  using 
smoothed  short  transmission  scan  data  are  shown  in  Fig.  5.23  (b),  5.23  (c),  and 
5.23  (d),  respectively.  From  this  figure,  we  can  observe  that  the  ML  method 
produces  an  image  with  more  details  than  the  other  images  (Note,  the  small  high 
density  spot  area  at  the  left  below). 

We  present  the  attenuation  images  of  the  ML  method  and  the  LS  reprojection 
method  in  Fig.  5.24  (b)  and  Fig.  5.24  (c),  respectively.  We  can  see  that  the 
ML  method  produces  a  more  accurate  attenuation  image  than  the  LS  reprojection 
method  fails  to  resolve  anything  (note  the  two  spots  at  the  center).  Compared  to 
the  attenuation  images  using  non-smoothed  transmission  scan  data,  attenuation 
images  here  have  more  details  and  features. 

The  comparisons  of  the  PSNR  for  the  emission  images  using  the  short  trans- 
mission scan  data  are  shown  in  Table  5.2.  Again,  the  ML  method  has  the  largest 
PSNR.  In  the  next  section,  we  show  the  results  using  the  real  data. 
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Figure  5.21.  Reconstruction  of  cardiac  attenuation  phantom  using  smoothed  short 
transmission  scan  data  (a)  attenuation  phantom  (b)  ML  method  (c)  LS  reprojec- 
tion  method 
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Figure  5.22.  ML  method  (a)  emission  phantom  (b)  emission  image  (c)  attenuation 
phantom  (d)  attenuation  image 
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Figure  5.23.  Reconstruction  of  chest  emission  phantom  using  smoothed  short 
transmission  scan  data  (a)  emission  phantom  (b)  EMML  algorithm  with  corrected 
data  (c)  LS  reprojection  method  (d)  ML  method 


84 


100 


100* 


v!™™-T-7r 


20 
40  i 
60  p 
80 
100^ 


20 


40  60  80  100 

(b) 


Figure  5.24.  Reconstruction  of  chest  attenuation  phantom  using  smoothed  short 
transmission  scan  data  (a)  attenuation  phantom  (b)  ML  method  (c)  LS  reprojec- 
tion  method 
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dB 

cardiac  phantom 

chest  phantom 

Algorithm 

Standard 

ML 

LS-Rep. 

Standard 

ML 

LS-Rep. 

Non-smooth 

-13.9 

22.13 

-9 

-5.8 

10.33 

9.28 

Smooth 

-13.62 

19.85 

-11 

-15.2 

9.02 

-8.22 

Table  5.2.  PSNR  (dB)  of  reconstruction  of  cardiac  emission  phantom  using  short 
transmission  scan  data 

5.2     Real  Data 

We  obtained  real  data  from  Dr.  John  Votaw,  Associate  Professor  of  Radiology 
at  the  Emory  University  Hospital.  To  create  the  data,  a  thorax  phantom  with 
a  hot  heart,  two  lungs,  a  spine,  and  soft  tissue  was  scanned  using  the  Siemens- 
CTI  ECAT  EXACT  921  PET  scanner.  This  scanner  is  comprised  of  24  rings  of 
diameter  82.5  cm,  where  each  ring  contains  384  detectors  of  size  6.35  x  6.35  mm2. 
In  addition,  the  diameter  of  the  scanner's  patient  port  is  56.2  cm  and  its  axial 
field  of  view  is  16.2  cm.  The  field  of  view  is  a  square  where  each  side  equals  to 
45  cm.  There  are  J  =  128  x  128  voxels  in  the  field  of  view.  In  two-dimensional 
mode,  virtual  rings  are  created  by  the  cross-ring  tubes.  Since  the  ECAT  EXACT 
scanner  consists  of  24  rings  and  23  virtual  rings,  it  provides  47  planes  of  data. 
Additional  details  on  the  scanner  can  be  found  in  the  paper  of  Wienhard  et.   al 

[53]. 

The  duration  of  the  emission  scan,  transmission  scan,  and  blank  scan  are  25 
minutes,  25  minutes,  and  60  minutes,  respectively.  Associated  with  each  plane 
are  the  sinograms  for  the  emission  data,  transmission  data,  and  the  accidental 
coincidence  data  from  the  emission  and  transmission  scan.  The  sinograms  are  of 
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size  160  x  192  corresponding  to  192  angles  with  160  bins  per  angle.  The  number  of 
tube  is  /  =  73536.  In  the  next  section,  we  discuss  details  of  the  implementation. 


5.2.1     Overview  of  Implementation 

All  the  images  shown  are  for  the  10th  plane.  This  plane  was  chosen  because 
it  contains  most  of  the  major  features  of  the  thorax  phantom.  The  10th  plane 
contains  0.288  million  emission  counts  (i.e.,  £(=1  j/,-  =  2.88  x  105)  and  13.82  million 
blank  scan  counts.  The  number  of  counts  for  the  low  quality  transmission  scan 
is  0.615  million  (i.e.,  £)Li  r<  =  6.15  x  105  ),  while  the  high  quality  transmission 
scan  contains  3.67  million  counts  (i.e.,  J2i=l  ri  =  3-67  x  1Q6)-  ^  snould  be  noted 
that  the  transmission  scans  were  automatically  smoothed  using  the  scanner's  data 
processing  algorithm.  The  emission  and  transmission  data  have  been  corrected 
for  detector  inefficiency  and  accidental  coincidences  using  standard  procedures, 
however,  no  scatter  correction  has  been  done. 

For  the  ML  method,  we  run  2  global  iterations,  7  and  5  iterations  for  the 
x-step  and  fi-step,  respectively.  In  other  words,  the  ML  method  requires  total  19 
iterations  and  24  iterations  to  get  an  emission  image  and  an  attenuation  image, 
respectively.  The  initial  estimate  for  .r-step  is  (x°)j  =  J2i=i  j%-  f°r  aU  h  where 
di  =  fi(^fc).  The  initial  estimate  for  /i-step  is  (fi°)j  =  0.07  for  all  j.  The  initial 
estimate  for  the  survival  probability  is  a°  =  gf^,  where  rb  =  60  minutes  and 
Tr  =  25  minutes. 

Twenty  iterations  of  the  EMML  algorithm  with  corrected  data  method  is  initial- 
ized using  x°-  =  ^-"f  v>  for  all  j,  where  yt  —  |S  For  the  LS  reprojection  method, 
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to  initialize  eq.  (4.1),  we  use  a  uniform  initial  estimate:  fij  =  0.07  for  all  j.  We 
run  20  iterations  to  get  he  LS  estimate  of  /*,  ft.  Then,  we  run  20  iterations  of 
the  algorithm  shown  in  eq.  (4.4).  Thus,  LS  reprojection  method  requires  total  40 
iterations  to  get  an  emission  image. 

The  transmission  scan  is  a  high  quality  scan  because  it  was  performed  long 
after  the  emission  scan  was  done  (no  emission  counts  will  be  mis-registered  as 
transmission  counts). 

5.2.2     Results 

The  emission  images  of  the  filtered  back-projection  method,  the  EMML  algo- 
rithm with  corrected  data,  the  LS  reprojection  method,  and  the  ML  method  are 
shown  in  Fig.  5.25  (a),  Fig.  5.25  (b),  Fig.  5.25  (c),  and  Fig.  5.25  (d),  respec- 
tively. From  this  figure,  we  observe  that  the  EMML  algorithm  with  corrected 
data  produces  an  image  with  clearer  lungs  boundaries  of  the  other  methods.  Also, 
the  ML  method  produces  comparable  but  smoother  image  than  the  EMML  algo- 
rithm with  corrected  data.  The  filtered  back-projection  method  provides  a  noisy 
emission  image  with  lots  of  artifacts  everywhere. 

Except  the  emission  images  of  the  filtered  back-projection  method,  all  emission 
images  have  a  "ring"  artifact  outside  the  phantom  (see  Fig.  5.25  (c)  and  Fig.  5.25 
(d)),  which  is  not  visible  in  the  image  using  the  EMML  algorithm  with  corrected 
data  because  the  limited  gray  level  scales  of  the  image  (the  heart  area  in  Fig.  5.25 
(b)  has  high  intensity  values  compared  to  the  background,  so  the  "ring"  artifacts 
on  the  background  cannot  be  seen). 
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Figure  5.25.  Emission  images  using  high  quality  smoothed  transmission  scan  data 
(a)  Filtered  back-projection  method  (b)  EMML  algorithm  with  corrected  data  (c) 
LS  reprojection  method  (d)  ML  method 


89 


The  attenuation  images  of  the  ML  method  and  the  LS  reprojection  method 
are  shown  in  Fig.  5.26  (a)  and  Fig.  5.26  (b),  respectively.  The  attenuation 
image  of  the  LS  method  in  Fig.  5.26  (b)  clearly  delineates  lungs,  spine,  and  tissue 
regions,  while  the  attenuation  image  of  the  ML  method  (see  Fig.  5.26  (a))  still 
suffers  from  line  artifacts.  Additionally,  the  attenuation  image  of  the  ML  method 
is  influenced  by  the  emission  data  (note  the  white  spot  at  the  center). 

We  also  show  the  emission  images  of  the  EMML  algorithm  with  corrected  data 
and  the  ML  method  using  the  non-smoothed  transmission  scan  data  in  Fig.  5.27 
(a)  and  Fig.  5.27  (b),  respectively.  The  result  of  EMML  algorithm  with  corrected 
data  becomes  very  noisy  while  the  image  of  the  ML  method  is  still  good. 

To  compare  the  computational  complexity,  we  can  observe  that  the  ML  method 
only  takes  19  iterations  to  obtain  its  best  image  while  the  LS  reprojection  and 
the  EMML  algorithm  with  corrected  data  requires  40  iterations  and  20  iterations 
to  obtain  the  best  emission  images,  respectively.  In  other  words,  the  ML  method 
produces  useful  emission  images  and  it  has  the  less  computational  complexity  than 
the  EMML  algorithm  with  corrected  data  and  the  LS  reprojection  method. 
In  the  next  chapter,  discuss  our  conclusion  and  possible  research  directions. 
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Figure  5.26.   Attenuation  images  using  high  quality  smoothed  transmission  scan 
data  (a)  ML  method  (b)  LS  reprojection  method 
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Figure  5.27.  Emission  image  using  high  quality  non-smoothed  transmission  scan 
data  (a)  EMML  algorithm  with  corrected  data  (b)  ML  method 


CHAPTER  6 
CONCLUSIONS 

In  this  dissertations,  a  ML  method  and  a  LS  reprojection  method  were  pro- 
posed for  correcting  attenuation  in  positron  emission  tomography.  The  ML  method 
simultaneously  estimates  the  emission  intensity  and  the  attenuation  density.  For 
the  LS  reprojection  method,  the  attenuation  density  is  first  estimated  using  a 
least  squares  estimator.  Then,  this  estimate  is  reprojected  to  estimate  the  sur- 
vival probabilities,  which  are  used  to  modify  the  probability  matrix.  Finally, 
the  emission  intensity  is  estimated  using  the  EMML  algorithm  with  the  modified 
probability  matrix. 

In  simulation  studies,  the  proposed  methods  were  compared  to  the  EMML 
algorithm  with  corrected  data  using  synthetic  data  and  real  data.  Below,  we  sum 
up  the  advantages  and  disadvantages  of  the  proposed  methods: 

•  The  ML  method  produces  emission  images  with  better  contrast  and  more 
details  with  less  computational  complexity  than  the  EMML  algorithm  with 
corrected  data  and  the  LS  reprojection  method  when  the  short  transmission 
scan  data  is  used. 

•  The  ML  method  provides  more  accurate  emission  images  than  the  EMML 
algorithm  with  corrected  data  and  the  LS  reprojection  method 
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•  The  LS  reprojection  method  produces  comparable  emission  images  but  is 
less  accurate  than  the  ML  method  while  the  transmission  scan  data  are 
smoothed. 

•  The  LS  reprojection  method  produces  better  attenuation  images  than  the 
ML  method  for  cardiac  attenuation  phantom. 

•  The  ML  method  produces  better  attenuation  images  than  the  LS  reprojec- 
tion method  for  chest  attenuation  phantom. 

•  The  EMML  algorithm  with  corrected  data  produces  much  better  emission 
images  when  the  transmission  scan  data  is  smoothed. 

Given  the  performance  of  the  proposed  methods,  it  is  felt  that  they  will  be 
most  useful  in  short  transmission  scan  protocols. 

6.1     Future  Work 

The  work  in  this  thesis  provides  a  new  approach  to  correct  the  attenuation  in 
PET.  With  the  simulations  using  the  real  data  obtained  from  the  Emory  Univer- 
sity, the  ML  method  steps  in  the  real  life  application  rather  than  the  theoretical 
work.  It  provides  an  alternate  tool  for  the  medical  society  to  use  the  research  result 
in  clinics.  With  more  research  done  on  searching  different  estimators  and  opti- 
mization algorithms  such  as  the  cyclic  coordinate  descent  golden  section  method, 
the  ML  method  can  be  expected  to  improve  the  emission  images  more.  In  ad- 
dition, we  may  use  the  information  from  the  histogram  of  the  reconstruction  of 
attenuation  density  to  segment  to  the  different  mediums.  Once  we  know  where 
the  different  mediums  distribute,  we  can  assign  the  common  attenuation  density 
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values  for  air,  water  and  bone.  Then,  we  can  use  that  information  to  get  a  better 
correction  for  attenuation. 

The  attenuation  density  reconstruction  can  be  used  to  locate  the  different 
tissues  such  as  bone,  water  and  air  within  ROI  by  setting  up  the  thresholds  on 
the  histogram  of  the  attenuation  density  reconstruction  images.  Then,  the  known 
accurate  attenuate  density  value  can  be  assigned  to  the  corresponding  (ij.  And 
fij  can  be  used  to  help  improve  the  emission  reconstruction. 

For  ML  method,  stopping  rules  for  the  rr-step  and  the  jz-step,  the  choice  of  the 
number  of  the  global  iteration  to  obtain  the  best  emission  images,  the  convergence 
of  the  algorithm,  and  the  statistical  analysis  of  the  estimate  are  open  issues.  For 
LS  reprojection  method,  a  stopping  rule,  and  the  analysis  of  the  algorithm  also 
needs  more  investigation. 


APPENDIX  A 
EMML  Algorithm  for  PET 

Shepp  and  Vardi  used  the  EM  algorithm  [16]  to  obtain  the  ML  estimate  of  the 
emission  intensity  for  error  free  data  [46].  In  their  formulation,  the  incomplete 
data  is  the  observed  data,  y,  and  the  complete  data  is  {N^},  where  Nij  is  the 
number  of  emissions  in  voxel  j  that  are  detected  by  tube  i.  Under  the  assumption 
that  {Nij}  are  independent  random  variables  and  the  assumptions  in  Section  2.2, 
the  EM  algorithm  consists  of  the  following  steps.  In  the  E-step,  the  expected 
value  of  the  log  likelihood  function  for  the  complete  data  given  the  incomplete 
data  and  the  current  estimate  of  the  emission  intensity  is  computed.  Let  L  be  the 
log  likelihood  of  the  complete  data: 

L(yn,Vi2,-  •  -,yu,  ■  ■  -,yn,yii,  ■  ■  -,yu)  =  log  P(Nn  =  yn,Ni2  =  yn,  (A.l) 
...,Nij  =  yu,...,Nn  =  yn,NI2  =  yi2,...,Nu  =  yu) .  (A.2) 

U(x,xk)     =     E[L(Nll,N12,...,N1j,...,Nn,NI2,...,NIj)\y,xk,x}     (A.3) 
/     J 
=      EiT,T,{-piJxJ  +  NijlogiPjXjWyiX^x]  (A.4) 

t'=i  i=\ 

=    E  £<-*;*;  +  fw$k*  h*W«»  •  (A-5) 

t=ij=i  \rmH 

In  the  M-step,  a  new  estimate  for  the  emission  intensity  is  obtained  by  maximizing 
U{x,xk)  with  respect  to  x: 

xk+1  =  arg  max    U(x,  xk)  .  (A.6) 
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Iterating  back  and  forth  between  the  E-step  and  M-step  leads  to  the  following 
iteration: 

where  xjf  is  the  kth  estimate  of  Xj.  This  iteration  is  called  the  expectation  max- 
imization maximum  likelihood  (EMML)  algorithm.  Details  on  the  EMML  algo- 
rithm can  be  found  in  [46]. 
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